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1.0 RESEARCH 


1.1 INTRODUCTION 

Probabilistic structural analysis methods are particularly useful in the design and 
analysis of space propulsion systems operating in severe and uncertain environments. By 
quantifying the uncertainties associated with the design, these methods play a critical role in 
establishing increased performance and reliability for such systems, both current and future. A 
prime example of a current space propulsion system is the Space Shuttle Main Engine (SSME). 
For this project, a fuel turbopump blade, a component of the SSME, was analyzed using a 
probabilistic finite element code, NESSUS, along with an embedded Material Strength 
Degradation (MSD) model. 

The use of NESSUS and the MSD model are dictated by the dispersion of the physical 
quantities involved. The considerable scatter of experimental data and the lack of an exact 
description of the underlying physical processes for the combined mechanisms of fatigue, 
creep, temperature variations and so on, make it natural, if not necessary to consider 
probabilistic models for a strength degradation model. This engine operates in a harsh 
environment with high loads and temperature extremes. As noted from C.C. Chamis and D.A. 
Hopkins “. . . hot rotating structural components are relatively small. Fabrication 
tolerances on these components which in essence are small thickness variations, can have 
significant effects on the component structural response. Fabrication tolerances by their 
very nature are statistical. Furthermore, the attachment of components in the structural 
system generally differ by some indeterminate degree, from that assumed for designing 
the component. In summary, all four fundamental aspects: (1) loading conditions, (2) 
material behavior, (3) geometrical configuration and (4) supports on which structural 
analysis are based, are of a random nature. The direct way to formally account for all 
these uncertain aspects is to develop probabilistic structural analysis methods where all 
participating variables are described by appropriate probabilistic density functions. [1 ]’’ 
The SSME fuel turbopump blade discussed here is a prime example of a small component 
whose variations in thickness, strengths, etc., significantly affect its structural response. 

The main objective of this project was to determine which input parameters most 
influence the structural response of the turbopump blade. A finite element model of a 
turbopump blade, composed of Inconel 718, a material for which the MSD model is calibrated, 
was analyzed for the effects of high temperature and high cycle mechanical fatigue. Section 
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1.2 presents NESSUS and the newly implemented MSD model. An overview of statistical 
terminology and definitions concerning structural reliability are provided in Section 1.3. 

Section 1.4 characterizes the SSME fuel turbopump blade as well as outlines NESSUS 
input and output data files. Sections 1.5 through 1.7 present analysis results. In Section 1.5 
reliability results are presented and comparisons between the Advanced Mean Value (AMV) 
method and the well-known Monte Carlo Simulation method are summarized. Section 1 .6 
discusses sensitivity factor results. Also discussed is the effect on sensitivity factor results for 
a change in the blade thickness random field from a single correlated random variable to a 
partially correlated random variable. Section 1.7 examines how changing the distribution type 
on a single random variable influences both reliability and sensitivity factor results. 
Conclusions drawn from these results are presented in Section 1.8. 
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1.2 NESSUS 

Historically, NESSUS has provided probabilistic structural analysis for the response of 
engineered structures. It integrates both probabilistic analysis and finite element methods to 
produce a useful tool for both design and analysis of critical structural components. NESSUS 
allows for both random and deterministic descriptions of input and output quantities. For 
example, input parameters for material properties such as Young’s modulus, Poisson’s ratio, 
and mass density can be modeled as random variables wherein the user provides the 
appropriate statistics, namely mean, standard deviation and distribution type [2]. In addition, 
random input can include random fields, as well as random variables. Distributions such as 
Gaussian, Weibull, or lognormal may be selected. Other input parameters, such as loading 
type, geometry, boundary conditions and initial conditions can also be considered random. 
Then, probabilistic analysis yields the cumulative distribution functions (CDF’s) of the output 
random variables, e.g., stresses, strains, displacements, etc. 

Code Organization 

NESSUS is divided into ten parts. The preprocessor, NESSUS/PRE, provides for the 
representation of random fields. NESSUS/PRE computes a set of independent random 
variables from the random fields using an eigenvalue decomposition [3]. PFEM is the module 
for coordinating the operations between the finite element code and FPI. Next, a deterministic 
finite-element module, NESSUS/FEM, generates the deterministic responses of the finite- 
element structural model to a number of prescribed perturbations of the input variables. 
NESSUS/FEM contains a mixed iterative procedure to obtain the response of a structure. A 
Newton-Raphson iterative procedure based on a mixed variational principle can be used in 
NESSUS/FEM. Besides FEM, BEM, the boundary element module is used for structural 
analysis. These deterministic solutions provide the data from which a first or second-order 
Taylor series of the response is fit that defines the explicit functional relationship between the 
input random variables and the response of interest. The postprocessor, NESSUS/FPI (Fast 
Probability Integration), performs the probability calculations. The postprocessor also includes 
a choice of probabilistic analysis methods other than FPI, including for example, Adaptive 
Importance Sampling (AIS) and Monte Carlo simulation [4,5]. SEMFEM however, contains 
the driving module for the Monte Carlo as well as the Latin Hypercube sampling of the 
structure. There are also other packages for specialized analysis such as the RISK, SYSTEM 
and SYSRSK modules. The RISK module computes the risk with respect to cost or a user 
defined criteria. The SYSTEM module is an alternate system reliability method for 
probabilistic fault tree analysis. The SYSRSK module computes the system risk by combining 
the cost of failure, in terms of replacement, inspection, repair or other criteria of individual 
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components of the system with their probability of failure. The current version of NESSUS 
provides these ten component parts in an integrated package. These modules can be used in 
various combinations as the user desires. 

This project required the use of the PRE, PFEM, FEM and FPI modules. For the 
random thickness field, PRE was used to compute a set of independent random variables from 
the random field thickness input using an eigenvalue decomposition. The resulting file was 
then used by the PFEM, FEM and FPI modules for further analysis. 

The AMV method is a mean-based predictor-corrector method which uses local 
gradients to project the results for the entire CDF then corrects the prediction with a re-analysis 
of the structural model at each predicted point. The local gradients are usually, but not 
necessarily, taken about the mean values of the random variables. The AMV method with 
iterations, termed AMV+, is an extension of the AMV method. It involves new gradient 
computations about the previously predicted results. These new gradients are then used to 
compute an improved result. Both AMV and AMV + results can either be specified as an 
automatic CDF analysis, a number of probability levels or a specified response corresponding 
to a specified probability [6], Convergence to a specific tolerance can be specified by the user. 
The FEM module then calculates the deterministic response where the Newton-Raphson 
iterative procedure is enabled. Finally the FPI module calculates the structural reliability of the 
turbopump blade. 

Material Strength Degradation fMSD) Model 

The enhanced version of NESSUS reported in this document, includes a previously 
developed probabilistic material strength degradation (MSD) model [7], The MSD model in the 
form of a postulated randomized multi-factor equation provides for quantification of uncertainty 
in the lifetime strength of components subjected to a number of diverse random effects. 
Presently, the model includes five effects that typically reduce lifetime strength: 

• high temperature 

• high-cycle mechanical fatigue 

• low-cycle mechanical fatigue 

• creep 

• thermal fatigue 

The MSD model was calibrated for INCONEL 718 by appropriate curve-fitted least squares 
linear regression of experimental data. Linear regression of the data for each effect resulted in 
estimates for the empirical material constants, as given by the slope of the linear fit. These 
estimates, together with ultimate and reference values, were used to calibrate the model 
specifically for Inconel 718 [8]. Lifetime material strength results, in the form of cumulative 
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distribution functions (CDF’s), illustrate the sensitivity of lifetime strength to any one, or a 
combination of, the five effects. The mathematical expression for this MSD model is: 



where, A it A i(J and A i0 are the current, ultimate and reference values, respectively, of a 
particular effect; is the value of the calibrated empirical material constant for the i ,h effect 

terms of the variables in the model; and S and S 0 are the current and reference values of 
material strength. Each term has the property that if the current value equals the ultimate value, 
the lifetime strength will be zero. Also, if the current value equals the reference value, the term 
equals one and strength is not affected by that value. The product form of equation (1.2.1) 
assumes independence between individual effects. This equation may be viewed as a solution 
to a separable partial differential equation in the variables with the further limitation or 
approximation that a single set of separation constants, a j5 can adequately model the material 
properties. 

The multifactor equation for material strength degradation (MSD) was implemented in 
NESSUS using the NZFUNC subroutine that was initially developed for predefined resistance 
models[9]. All input from the MSD model is in the NESSUS/PFEM input deck using a 
keyword interface consistent with previous versions of NESSUS. 

The model is currently set up for five fixed effects and twelve general effects that can be 
defined by the user. The fixed effects include: Temperature, High cycle mechanical fatigue. 
Low cycle mechanical fatigue, Creep, Thermal fatigue. When expanded for all effects, the 
equation is as follows: 



Temperature HCF LCF Creep Thermal Fatigue General 


where T is temperature, N is cycles, t is time, and A is a user-defined effect. The lower case 
exponents are the empirical material constants for each effect. The subscript, u, refers to an 
ultimate value, the subscript, o, is the reference value, and the non-subscripted term is the 
current value for the effect. 

To increase model sensitivity for the effect, a logarithmic base ten transformation is 
introduced. For Inconel 718, all effects except temperature require the logarithmic 
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transformation and have the following form. 


S = 


LOGjA^-LOGW 
LOG{A iu ) - L0G(A jO ) 


(1.2.3) 


For other materials, the nature of the data will dictate the use of the log transformation. 
The effects to be used and the model type (log transformation) are defined in the g-function 
definition section of the PFEM input file. Any combination of effects can be selected and any 
of the terms can be considered random if desired. The form of the g-function used by 
NESSUS is 




“ < 7 , 


(1.2.4) 


where S 0 is the reference value of material strength, A i(J is the ultimate value of the particular 
effect, A l0 is the reference value of the particular effect. A, is the current value of the particular 

effect, a, is the empirical material constant for the particular effect, and a is the structural 

response as calculated by NESSUS/FEM, i.e. stress. The equation expanded for the two 
effects of high temperature and high-cycle mechanical fatigue is shown in equation 1.2.5 
below: 


z 


l r.-rT 

T LOG(Ny)- LOG(N t ) " 

» 

1 

fx* 

i 

LOG(N u )-LOG{N 0 ) 


(1.2.5) 


With the implementation of the MSD model, NESSUS can now determine the reliability of a 
structural component utilizing its material strength variability. Thus for example, NESSUS 
now produces statistical distributions for the component material strength, S, as well as a 

statistical distribution of component stress, a. 
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1.3 RELIABILITY 


Structural reliability calculations can be made using the Material Strength Degradation 
(MSD) model combined with the stress output from the FEA model. Structural reliability is the 
probability that the strength exceeds the stress as defined by equation 1.3.1 below: 

R = P(S><7 ) , (1.3.1) 


where R is the structural reliability and P(S > o) is the probability that the strength is greater 

than the stress. The von Mises stress term is used in this project since it represents a three- 
dimensional state of stress and is considered a good indicator for the onset of mechanical 
failure for ductile materials [10]. The equation for Von Mises stress is given below: 

(tr, -ct 2 ) 2 +(ct 2 -ct 3 ) 2 +(<r, -cr 3 ) 2 


<7 = 


-. 1/2 


(1.3.2) 


where a,, a 2 and a 3 are the principal stresses. These principal stresses are the roots to 

the following cubic equation [11]: 

o' 3 " K + ^ y + (<7,cr v + a x o z + <r v a z -t 2 „- < - )tr 

' (1.3.3) 

In short, von Mises stress is a scalar value representing a three-dimensional state of stress. 

The distortion-energy theory predicts that yielding will occur whenever the distortion energy 
per unit volume equals the distortion energy in the same volume when yielded in a simple 
tensile test as shown in equation 1.3.4 below [12]. 

o' > S Y (1.3.4) 

The reliability approach is a departure from the factor of safety approach widely used in 
engineering design. Input of the usual physical quantities as either random variables or random 
fields has many advantages including: 

• probability of failure- quantifies reliability 

• sensitivity-identifies important variables 

• importance 

Using the probabilistic approach, an analysis of the turbopump blade model yielded a 
reliability of 0.999. Initially all of the random variables were modeled as normal distributions 
and the blade thickness was modeled as a fully correlated random field. Taking the median 
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value of the strength and the deterministic value of the stress from NESSUS, the factor of 
safety is: 


N = 


5 ; 

o 


1.16107E05 
1.0554 £05 


(1.3.5) 


The reliability measure of 0.999 denotes a failure for 1 in every 1000 uses. The factor 
of safety, N=1.10, implies the structure can withstand 1.1 times its nominal load. This may 
give the engineer a false sense of security. The factor of safety approach requires engineering 
judgement that results from experience. 

The structure is also under high cycle fatigue. To this end the Goodman Equation is 
used [13] as given by Equation 1.3.6: 


£- 



1 

n 


(1.3.6) 


where a m is the mean stress due to cyclic loading, o a is the average stress, S U1 is the ultimate 

strength of the material and S e is the endurance limit of the material. This results in a factor of 
safety of 1 . 1 1 as shown below. 


32.5 ksi 105.5 ksi T 1 

111 ksi 173 ksi 


(1.3.7) 


Reliability has advantages in that it quantifies the uncertainty of material properties and 
provides for the driving variables for design. While the reliability method has its advantages, it 
is at the expense of being more complex requiring a background in statistics and the 
comprehension of terms such as the performance function. 

The performance function is a user-defined function. In this particular case, the 
performance function is defined as the equation: 


Z = S - a, (1.3.8) 

where Z is a measure of performance such as stress, strain or displacement, S is the strength 
and a is the von Mises stress [14]. The strength value resulting from the combined effects of 

high cycle fatigue and high temperature was defined previously in the MSD model. The von 
Mises stress, is used to compare with material strength and estimate failure or nonfailure. 

The limit state is the surface that separates the density space into safe and fail regions 
where the performance function is zero. This limit state is termed g, as in the 
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following equation: 

g(JO = SW-<T(X) = 0 . (1-3.9) 

where g(X) is the limit state, defined on the domain of the vector of random variables X, S(X) 
is the strength and a(X) is the stress. When a and S are independent, i.e. no common random 

variables, their two distributions are as shown in Figure 1.3.1. Note that there is some overlap 
between the stress and the strength. The region where this overlap occurs is the probability of 
failure. For most design cases, attempts are made to minimize this overlap for a smaller 
probability of failure. 

Within the vector X of equation 1 .3.9, there are 13 random variables and 1 random 
field. These are further divided into 9 random variables for the MSD model and 5 for the stress 
side. The MSD model utilized includes the effects of: High-Cycle Fatigue (HCF) and High 
Temperature. The stress side’s random variables are the loading, modulus of elasticity, 
Poisson’s ratio, mass density and the thickness random field. The HCF components are the 
ultimate number of cycles, designated number of cycles, reference number of cycles and an 
empirical material constant (exponent). The temperature side of the strength includes the 
ultimate temperature, designated temperature, reference and an exponent. The remaining 
variable is the material’s initial strength. These 13 random variables and 1 random field result 
in a complex failure hypersurface which is not easy to visualize. Fortunately NESSUS 
resolves this problem for the user. 
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For the case of only two random variables included in a structural reliability analysis, 
the Joint Probability Density Function (JPDF) can be visualized. The probabilities are contours 
that are arranged such that the lowest probability is the biggest ellipse and the highest 
probability is the smallest ellipse. Thus the highest probability occurs at the center of the 
innermost ellipse in Figure 1 .3.2. Recall that the limit state function is where the performance 
function equals zero, i.e. g = 0. This limit state function is usually desired to be some distance 
away from the location of the highest probability as shown in Figure 1 .3.2. The probability 
that the structure will fail is upward and to the right and above the performance function. 
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Structural reliability is the probability that the structure will not fail. This is given by 
the following equation: 

Reliability = P(g > 0) (1.3.10) 

where P is the probability. NESSUS output results yield the probability of failure, p f . In 

terms of this failure, the reliability is calculated from equation 1.3.1 1 below: 

Reliability = 1 - P(g < 0) , (1.3.11) 

where P(g<0) is the probability of failure. Thus, the area to the right of the origin in Figure 
1 .3.3 is the reliability of the particular structure and the shaded area to the left is the failure 
region. 


PDF 



Figure 1.3.3 Overlap of Strength - Stress 
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Most Probable Point (MPP) 

The Most Probable Point (MPP) is the point on the limit state with the highest joint 
probability as shown in Figure 1 .3.4. Note that this joint PDF is three-dimensional and is now 
centered about the origin. This is due to the Rosenblatt transformation. The Rosenblatt 
transformation can be used for correlated random variables and is a transformation whereby 
distributions are transformed from one distribution to another [ 15 ]. In NESSUS, the 
Rosenblatt transformation is impractical because the available data is often insufficient to 
establish the joint and the conditional probability distributions In this case the usual random 
variables are transformed from an original input distribution into a standard normal 
distribution. This is where the mean of the distribution is zero and the standard deviation is 
one. 



Cut volume = P(Z<=0) 


Figure 1 .3.4 Limit State in u-Space for Two Random Variables 
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Beta 


As the random variables are transformed into the standard normal distribution, the 
resulting distance from the origin to the MPP, termed beta, is defined by equation 1 .3. 12 . 


P 2 = K+K+K +•••(«)! 

= 2 >>? 


(1.3.12) 


i=i 


where the u, are the transformed random variables in standard normal space at the limit state. 
The probability of failure is approximately the value of the joint probability at the MPP as seen 
in the equation below: 


P, = <K/S) = J..f ,JMdu (1.3.13) 

where / u is the joint probability distribution function and g approx is the approximation of the 
limit state g. This approximation of g can be a polynomial, usually either a linear or quadratic 
equation. 


Sensitivity Factors 

The sensitivity factors are a ranking of which random variables “drive” or most 
influence the output results, e.g. probability of failure. These are useful in the assessment of 
which random variables need to be changed for design optimization. 

These sensitivity factors are defined as the direction cosines to the MPP of the input 

random variables for NESSUS. The sensitivity factor, a, is given by equation 1.3.14 below. 


a.\ 


fdg} 
V dX i J 


(1.3.14) 


dg 


where a, is the sensitivity factor, i indexes the random variable number, is the gradient of 
1 ' OX; 


the limit state with respect to the particular random variable Xj, and a, is the standard deviation 


of the ith random variable. When the limit state is taken as a linear approximation, the 
following equation holds: 


a. = 


■fowl 

■?(*,) 


a 






(1.3.15) 


where g is the linear approximation of the limit state, n is the number of random variables, and 
a c and a ; are constants determined by Taylor’s series expansion NESSUS [16]. 
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Distribution Types 

The lognormal and Weibull distributions are derived from the normal distribution. The 
lognormal distribution was calculated using the following equations [17]: 

C x = — . ( 1 . 3 . 16 ) 


where C x is the coefficient of variation, G x is the standard deviation and |i x is the mean from the 
normal distribution. Likewise for the median the equation is: 


X = 




i/i+cf ' 

where is the median of the distribution. 

The Weibull distribution was calculated using the following equations: 


( 1 . 3 . 17 ) 


f x = 


where a is defined as; 



a-1 

" 


a " 

a 

X 


X 


kp) 

p exp 


Ifi) 



x>0, a >0, /? > 0 , 


( 1 . 3 . 18 ) 


P is defined as; 



and T is defined as the integral; 


T(jc) = Jr'e-'dt 

o 


( 1 . 3 . 19 ) 


( 1 . 3 . 20 ) 


( 1 . 3 . 21 ) 


Monte Carlo Simulation 

Reliability results were calculated using the Advanced Mean Value Method (AMV) and 
the well-known Monte Carlo Simulation method. Using the Advanced Mean Value Method 
requires several steps: computing sensitivities by perturbing each random variable and 
recomputing the response; using the Fast Probability Integration (FPI) to compute the Mean 
Value (MV) cumulative distribution function (CDF.) and associated Most Probable Point 
(MPP); the response at the predicted MPP for each probability level; and computing the 
sensitivities about the predicted MPP. The last two steps are AMV and AMV+ respectively. 
Once these last two steps are performed, the run proceeds through FPI to compute the response 
at each probability level and new MPP [18]. With the FPI calculated the convergence criteria is 
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measured at each probability level. If the response does not converge, the data is run through 
the AMV and AMV+ steps until convergence is reached [19]. 

Monte Carlo Simulation is a more widely accepted method for calculating reliability 
results. The method is more robust than AMV+ in that it is independent of the number of 
random variables. Monte Carlo Simulation uses random samples from the probability 
distributions of each random variable and computes the response. The probability of failure is 
determined by counting the number of failures and dividing by the total number of samples. 
The method is exact as the number of simulations approaches infinity. The major disadvantage 
of this method is the computing time required for large sample sizes. Thus the number of 
simulations (samples) must be given careful consideration. In many “quick and dirty” 
calculations, a very simple equation for determining the sample size is shown in equation 
1.3.22 below [20]: 

N > — (1.3.22) 

Pf 

where, N is the number of samples for the Monte Carlo simulation required for the given 
probability of failure; p f The Shooman equation gives a more accurate assessment of how 
many points are necessary for a given probability of failure. The Shooman equation is that 
given by Equation 1.3.23 [20]: 


%error = 100 • O '[(1 - a)/2] ■ 



(1.3.23) 


where (1- a) is the confidence desired and 0~'[(l-a)/2] is the inverse probability of the 

confidence divided by two. For calculating the number of points for a 95% confidence the 
Shooman equation is as follows: 



Pf 


200 

%error 


(1.3.24) 


where the %error could be some reasonable value, usually between 5% and 20%. 


NASA/CR— 2001-211112 


15 




1.4 TURBOPUMP BLADE AND DATA FILES 

Previously, the probabilistic analysis of a single structural member, a MAR-M246 
SSME turbopump blade was undertaken [21]. The finite element model used is a reasonably 
realistic representation of the Space Shuttle Main Engine fuel tuibopump blade shown below in 
Figure 1.4.1. 


Platform 


Trailing edge 


Firtree or root 



Airfoil pressure side 


Leading edge 


Figure 1.4.1 Drawing of the Physical Turbopump Blade 


Note the overall shape of the structure that was modeled. It has a curved shape toward 
the tip of the aerodynamic surface and a less noticeable twist along the length of the turbopump 
blade. The various parts of the blade should be noted. They are the firtree root, the platform, 
the trailing edge, the leading edge, the airfoil suction side and the airfoil pressure side. The 
firtree root, sometimes called just the root attaches the turbopump blade to the rotor. Above 
the root is the platform. It separates the root and nonaerodynamic parts of the rotor from the 
fluid flow of the blade. Important to the fluid flow are the leading edge and the trailing edge. 
These parts are where the fluid flows toward and away from the blade, respectively. As 
mentioned before, this blade has curvature like a wing. Hence, to one side of the blade there is 
an area of high pressure and low velocity relative to an area with low pressure and high 
velocity. The high pressure side is called the airfoil pressure side and the low pressure side is 
called the airfoil suction side. At the right side, note that the thickness of the blade starts to 
decrease from the leading edge to the trailing edge. This is an important attribute to keep in 
mind when analyzing the resulting stresses. Figure 1 .4.2 shows the boundary conditions and 
forces on the blade. 


N AS A/CR— 2001-211112 


17 



Note that on the left-hand side of the blade there is a cantilever on the center node 
surrounded by trolleys on the four outer nodes. These outer nodes are constrained to move 
transversely in the Z-direction so that the blade is allowed to twist within its attachment on the 
rotor as the torque of the aerodynamic force loosens the roots fit. Other forces such as 
pressure and thermal gradients are neglected. This blade is under a centrifugal load due to a 
rotational speed of 37,000 rpm according to Milton L. Littlefield of Pratt and Witney, “Mission 
conditions that include shaft speeds of 37,000 rpm, temperature extremes ranging from minus 
265 to 2,300 °F and pressures over 8,000 psi place high stresses on the hardware”[22]. 

The finite element model of the SSME turbopump blade consists of 55 nodes and 40 
four-noded isoparametric elements. Since this is a relatively thin blade, the elements are shell 
elements. These elements consist of 4 nodes and their corresponding nodal thicknesses 
composed of 5 layers. 



Figure 1 .4.2 Drawing of the Blades Constraints and Forces 
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Random Variable Inputs for Stress 

Table 1.4.1 lists the Inconel 718 random variable inputs for the stress calculations. 
Inconel 718 is a precipitation-hardenable nickel-chromium alloy containing significant amounts 
of iron, niobium and molybdenum along with lesser amounts of aluminum and titanium. 
Inconel 718 has excellent creep-rupture strength with properties of high strength, corrosion- 
resistant material suitable for temperatures from -423 to 1 300 °F [23]. 

These random variables include the modulus of elasticity [24], Poisson’s ratio, mass 
density and centrifugal loading . The random variables above have defined means and standard 
deviations. The coefficients of variation are initially assumed to be 5% as shown by the 
equation 1.4.1. Initially, the distribution type for every random variable was assumed to be 
normal. 


Std.dev. _ 1.240£03 
mean 2.480£04 


(1.4.1) 


Table 1.4.1 Random Inputs for Stress 


FINITE ELEMENT RANDOM VARIABLES 

RANDOM 

VARIABLE 

NAME 

MEAN 

STANDARD 

DEVIATION 

DISTRIBUTION 

MODULUS OF 
ELASTICITY 

EMOD 

2.48E04 (ksi) 

1.24E03 (ksi) 

NORMAL 

POISSON’S RATIO 

PMOD 

0.271 

0.01355 

NORMAL 

MASS DENSITY 

DMOD 

0.296 (lbm)/in 3 

0.0148 (lbm)/in 3 

NORMAL 

CENTRIFUGAL 

LOADING 

DLOAD 

37,000 (rpm) 

1850 (rpm) 

NORMAL 

BLADE THICKNESS 

THICK 

0.24" - 0.097" 

0.012" - 0.0005" 

NORMAL 
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Correlated Random Fields 

The nodal blade thicknesses were input as a random field. Initially this was done 
initially to compare to the previously run problem [21]. Two different random field analysis 
were conducted: a fully correlated random field and a partially correlated (95% uncertainty 
retained). The partially correlated random fields are those where the variance-covariance matrix 
will be densely populated but diagonally dominant. The correlation coefficient is defined by 
the following exponential decay equation: 

Po = * ' . (1.4.2) 

where p tJ is the correlation coefficient, l is the characteristic length, X ( and X, are the 

coordinates of point i and j respectively. The characteristic length is estimated measured or 
related to other physical (measurable) characteristics. This makes for the following variance- 
covariance matrix: 

^x,x 1 = Pifix^Xj ’ (1.4.3) 

where p jJ is the correlation coefficient, a Xj and a Xj are the standard deviations of the field 

variable X at points i and j respectively. One item to note is that the transformed random 
variables are the eigenvectors of the variance-covariance matrix. The associated eigenvalues 
may be interpreted as the square of the standard deviation of the transformed variables. 

For this model, the characteristic length is 1 inch. Since the length of the blade is 1 .25 
inches and a mesh size is 0. 1 inch by 0.2 inch, the characteristic length of 1 inch is plausible. 
The uncertainty retention measure determines how many quantities need to be retained to 
capture the randomness of the structure and to shorten the computational time. This means that 
for small standard deviations, these uncertainties may be neglected to substantially decrease the 
number of random variables. For certain strongly correlated random fields it is not uncommon 
that only a few of the significant random variables account for the majority of the uncertainties 
in the random field. Inside of NESSUS this retained uncertainty is an input to NESSUS/PRE. 
The fully correlated random fields are those where the variance-covariance matrix is 

singular and fully populated. The correlation coefficient p 0 is equal to 1, thus having the 
entries of the variance-covariance matrix equal to the equation: 

Cx/Xj = G x o x . , (1.4.4) 

where cr Xi and cr Xj are the standard deviations of the field variable X at points i and j 
respectively [25]. 
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Random Variable Inputs for Strength 

The random variable inputs for strength are listed in Table 1.4.2. Note that these 
random variables consist of three different types: High-Cycle Fatigue (HCF), High 
Temperature and Strength. As stated previously, the random variables for HCF include 
ultimate number of cycles, designated number of cycles, reference number of cycles and the 
empirical material constant. Likewise, variables for the effect of high temperature include the 
ultimate temperature, current temperature, reference temperature and empirical material constant 
for temperature. The last random variable is the initial strength. Note that like the random 
variables for stress, the standard deviation is 5% of the mean and the distribution type was 
initially assumed to be normal. 

A sample input file for NESSUS is given in the Table 1.4.3. It starts with a keyword 
called *PFEM which stands for Probabilistic Finite Element Method. The next major keyword 
is the *ZDEFINE which defines the definition of the performance function. 
*EXPLICITMODEL defines the input of random variables for the MSD model. The keyword 
*ZFUNCT 3 0 defines which reliability model to use, with 3 corresponding to the MSD 
model. Below this is the *MSDM input, which is for the Material Strength Degradation model. 
The particular effects and their input are placed here. The linear or log type is also designated 
in the MSDM section. The keyword *COMPUTATIONALMETHOD designates the number 
of random variables on the stress side. They are given a random variable number and name 
designation below the *RVDEFINE keyword. The COOR is the thickness change that is 
applied to the coordinates. The first number is an integer defining the nodal number. 

Following that is the x, y, z and thickness coordinates. Note that the x, y, and z coordinates 
are zero. This means that the only change applied to the coordinates is a change in thickness. 
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Table 1.4.2 Random Inputs for Strength 



PROBAE 

H Li STIC RANDOM VARIABLES 

TYPE 

RANDOM 

VARIABLE 

SYMBOL 

MEAN 

STANDARD 

DEVIATION 

DISTRIBUTION 

TYPE 

HCF 

ULTIMATE 
NUMBER OF 
CYCLES 

NU 

1.0EI0 

5.0E08 

NORMAL 

HCF 

DESIGNATED 
NUMBER OF 
CYCLES 

N 

2.5E05 

1.25E04 

NORMAL 

HCF 

REFERENCE 
NUMBER OF 
CYCLES 

NO 

0.25 

0.0125 

NORMAL 

HCF 

EMPIRICAL 

MATERIAL 

CONSTANT 

W 

0.141 

0.00705 

NORMAL 

High 

Temperature 

ULTIMATE 

TEMPERATURE 

VALUE 

TU 

2369 °F 

118.45 

NORMAL 

High 

Temperatire 

CURRENT 

TEMPERATURE 

T 

i j 

1000 »F 

50.0 

NORMAL 

Higi 

Temperature 

REFERENCE 

TEMPERATURE 

TO 

75 oF 

3.75 

NORMAL 

High 

Temperatire 

EMPIRICAL 

MATERIAL 

CONSTANT 

QQ 

0.2422 

0.01211 

NORMAL 

Strength 

INITIAL 

STRENGTH 

SO(psi) 

I.480E5 

7.40E03 

NORMAL 
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Table 1.4.3 Sample NESSUS Input File 


C ssmehcf718nrvla.dat The input file for NESSUS 

* PFEM 
C 

C Z- FUNCTION 

C 

* ZFDEFINE 

* EXPLICITMODEL 9 

6, 7 , 8, 9, 10, 11, 12, 13, 14 

* ZFUNCT 3 0 
* MSDM SO 

C Definitions of input for the Random 

TEMP LINEAR TU T TO QQ C Variables of the Material Strength 

HCF LOG NU N NO W C Degradation model. 

END 

* COMPUTATIONALMETHOD 15 C Random Variable memory allocation for the 

12345 C stress side. 

* END 

C RANDOM VARIABLE DEFINITIONS C Defines the random variables. 

* REDEFINE 

* DEFINE 1 C Defines the first random variable, a random field. 

C The random field is run through the NESSUS 
C preprocessor , NESSUS/PRE. 

C While this is in the form of a coordinate input deck, it has zeros in 
C the x, y and z coordinates. The only nonzero entry is the thickness. 

C This is the change that is applied to the coordinates in the blade. 
THICK1 

1.3650E-01 6.8249E-03 NORMAL 

COOR 

1 0 . 0000E+00 0 . 0000E+00 O.OOOOE+OO 1.7583E+00 

2 0.0000E+00 O.OOOOE+OO O.OOOOE+OO 1.7583E+00 

3 0.0000E+00 O.OOOOE+OO O.OOOOE+OO 1.7583E+00 

4 0.0000E+00 O.OOOOE+OO O.OOOOE+OO 1.7583E+00 

5 0.0000E+00 0.0000E+00 O.OOOOE+OO 1.7583E+00 

6 O.OOOOE+OO O.OOOOE+OO O.OOOOE+OO 8.0001E-01 

7 0.0000E+00 0.0000E+00 O.OOOOE+OO 1.4446E+00 

8 O.OOOOE+OO O.OOOOE+OO 0.0000E+00 1.6685E+00 

9 0.0000E+00 O.OOOOE+OO O.OOOOE+OO 1.7490E+00 

* DEFINE 2 C The second random variable. 2 is its random number. 

DMO D 

2.48E07 1.24E06 NORMAL 

PROP 75 

1 40 0 1 0 0 0 

* DEFINE 3 

DMOD 

0.271 0.01355 NORMAL 

PROP 75 

1 40 0 0 1 0 0 

* DEFINE 4 

DMOD 
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Table 1.4.4 contains an excerpt from the .mov output file. This file includes a 
summary of the results from the probabilistic analysis of the blade. The first line in Table 
1.4.4 shows probability of failure as 0.00022074 corresponding to the limit state of Z=0 at 
-3.5140 standard deviations from the mean. Below this line is a list of the random variable 
names, their values at the MPP, their corresponding deviations from their means and their 
sensitivity factors. Finally, the bottom line is a reprint from the top line with additional decimal 
places showing. 


Table 1.4.4 Sample NESSUS Output File 


C Filename: ssmehcf 7 1 8nrvla . mov A Sample output file from NESSUS. 

Response/Probability level: 1 

Level Z-value u(std. normal) Probability 

1 0 . 0 00 000E+00 -3.51401 0.000220740 

Most probable point (MPP) or design point 

Level 1: (Z-value = 0.00000E+00, u - -3.5140, Probability = 

0.000220740) 


R.V. name 

THICK1 

EMOD 

PMOD 

DM0 D 

DLOAD 

NU 

N 

NO 

W 

SO 

TU 

T 

TO 

QQ 

Z 

0 . OOOOOOOOE+OO 


X-value Std. 

0 . 1 309896E+00 
0 . 2 4 7 908 4E + 08 
0 . 2706960E+00 
0 . 8082256E-03 
0 . 4402004E+04 
0 . 99 9414 9E+1 0 
0 . 2502581E+06 
0 . 2498878E+00 
0. 1 422 7 99E + 00 
0. 1355315E+06 
0 . 2335924E+04 
0 . 101 4 30 1E + 04 
0 . 7495304E+02 
0 .2446292E+00 

CDF RESULTS 
U 

-0. 35140116E+01 


Dev. from Mean 
-0 . 807400 
-0.007384 
-0.022434 
1 . 101191 
2 . 722211 
-0.011701 
0.020652 
-0.008976 
0.181542 
-1.684931 
-0.279243 
0.286011 
-0.012522 
0.200592 


PROBABILITY 
0 . 22074004E-03 


Sensitivity factor 
-0 . 229766 
- 0 . 002101 
-0.006384 
0.313371 
0.774673 
-0.003330 
0 . 005877 
-0.002554 
0 . 051662 
-0.479489 
-0.079465 
0 . 081392 
-0.003564 
0.057084 


ITER. NO. 
8 


The mov output from NESSUS contains several important results. There is a title line 
followed by a description of what kind of output. The Level term is the level as defined by the 
user in the FPI input or is a predefined number of levels that NESSUS automatically chooses 
for a full CDF analysis. The Z-value is the performance level, followed by its location in 
standard normal coordinates and finally the probability of having that performance level. 
Following that is the MPP and another line containing the Z-value, the standard normal value 
and the probability. Next there is a line defining the random variable name, the x-value or the 
value of that particular random variable in normal space as given in the input file, the standard 
deviation of that x-value from it mean and finally its sensitivity factor. Finally, the last line 
below all the random variables is the same as the top of the input file except for the inclusion of 
more decimal places. 
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1.5 RELIABILITY RESULTS 


AMV+ and Monte Carlo Simu lation Results 

Table 1.5.1 compares the reliability results from the AMV + technique with the well 
known Monte Carlo Simulation method. Percent differences between the two methods were 
calculated. The number of samples for the Monte Carlo Simulation were increased from 
50,000 and 200,000 sample points. The run times varied from one day of CPU time for 
50,000 sample points, to six days for 200,000 sample points. 

Case 1 was for a fully correlated random field and 50,000 sample points. The 
probability of failure calculated for Monte Carlo was 2.60E-04 while the probability of failure 
from the AMV+ method was 2.207E-04. The percent difference was calculated to be 1 5. 1 %. 
Increasing the number of samples to 200,000 as shown in Case 3 reduced the difference to 
7.7%. This is in agreement with the reduction in percent error given by the Shooman equation 
from approximately 60% to 30% error, roughly one half. 

The additional number of samples increased the Monte Carlo run time from one day to 
six days. The AMV+ solution ran in 3 to 4 minutes, thus giving a substantial savings in 
processor time. With almost instant feedback from the AMV + and accuracy comparable to that 
of Monte Carlo, the AMV+ method offers the ability to examine more variables in a given 


amount of time. 

Case 2 was for the blade thickness modeled as a partially correlated random field and 
50,000 sample points. Monte Carlo resulted in a probability of failure of 3.40E-04, while 


AMV+ yielded a probability of failure of 2.207E-04. The percent difference between these two 
values was 33.2.%. As before, this difference was decreased by increasing the number of 
sample points used in the Monte Carlo Simulation from 50,000 to 200,000 sample points as 
seen in Case 4. The percent difference decreased from 33.2% to 14.7%. Thus, the AMV+ 
method yielded comparable results with much greater efficiencies than the Monte Carlo 
method. 


Table 1.5.1 Comparison of Probabilities of Failure and Reliability 


Case 

Number 

Blade thickness 
uncertainty 

Number of 
Samples 
for Monte 
Carlo 

Probability of 
Failure 

(Monte Carlo) 

Reliability 
Based on 
MC 

Probability 
of Failure 
(AMV+) 

Reliability 
Based on 
AMV+ 

%Error 

(Shooman) 

Companson 
%Difference 
AMV+ vs. MC 

1 

1 fully correlated 
random field 

50,000 

2.60E-04 

0.99974 

2.207E-04 

0.99978 

60.19 

15.1% 

2 

95% unoertanity 
retained .partially 
correlated 

50,000 

3.40E-04 

0.99966 

4.530E-04 

0.99955 

42.01 

33.2% 

3 

1 fully correlated 
random field 

200,000 

2.Q5E-04 

0.99980 

2.207E-04 

0.99978 

30.10 

7.7% 

4 

95% unoertanity 
retained .partially 
correlated 

200,000 

3.95E-04 

0.99961 

0.00045 

0.99955 

22.31 

14.7% 
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Reliability Results Overlap 

Figure 1 .5. 1 is a plot of the probability density functions (PDF’s) for both the resulting 
strength and stress distributions. First only the MSD model was run and then only the stress 
side of the SSME fuel turbopump turbopump blade was run. Next, the numerical derivatives 
of the MSD side were calculated and plotted. Finally the numerical derivatives of the stress side 
were calculated and plotted. The derivatives are given by equations 1 .5. 1 and 1 .5.2, where the 

p’s are the cumulative probabilities, the S’s are the strengths and the o’s are the stresses. 


dp _ P^zPq 
dS S l -S a ’ 

dp _ 

da cr, -a 0 


(1.5.1) 

(1.5.2) 


The overlap of these two PDF’s in Figure 1 .5. 1 is the region where the probability of failure 
exists. 


Overlap of Stress and Strength Plots for 

Inconel 718 @1000 °F for Blade Thickness 
as a Fully Correlated Random Field 



Figure 1.5.1 Overlap of the Strength and Stress Plots 
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Reliability results for Inconel 7 1 8 at a temperature of 1000 °F are shown in 
Figure 1.5.2. It shows the cumulative distribution functions CDF’s for the cases of fully 
correlated and partially correlated (95% uncertainty retained) random fields. This overlay on 
normal probability paper shows that there is very little difference between the probability 
results obtained for a fixlly correlated random field and those for a partially correlated random 
field. In the region where the blade fails, i.e. where (strength - stress) = 0, the shift is very 
minor. The probabilities for the fully correlated and the partially correlated (95% uncertainty 
retained) fields are essentially right on top of each other. 


Reliability for Inconel 718 @ 1000 °F 
Single Fully correlated and 95% Uncertainty Retained, 
Partially correlated Random Reids 



Figure 1.5.2 Reliability Plots on Normal Probability Paper 
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1.6 SENSITIVITY FACTOR RESULTS 

Table 1 .6. 1 states the sensitivity factors at failure for the SSME fuelpump tuxbopump 
blade with the thickness modeled as both a fully correlated random field and a partially 
correlated random field with 95% of the uncertainty retained. 

Note that the order ranking is for the fully correlated random field. For both fully and partially 
correlated random fields the centrifugal loading (DLOAD) has the highest sensitivity factor, 
followed by the initial strength (SO). Mass density (DMOD) and thickness (THICK) are third 
and fourth depending on correlation. The next four random variables, current temperature (T), 
the ultimate temperature (TU), the temperature empirical material constant (QQ) and the HCF 
empirical material constant (W), have sensitivity factors that are an order of magnitude smaller 
and are considered insignificant in comparison with other variables. The last six random 
variables are smaller by two orders of magnitude and are clearly insignificant in influencing the 
output. 


Table 1.6.1 Ordered Sensitivity Factors 


Random Variables 

Name 

Sensitivity 
Factor, a 
Fully Correlated 

Partially 
Correlated 95% 
Uncertainly 
Retained 

Partially 
Correlated 80% 
Uncertainity 
Retained 

centrifugal loading 

DLOAD 

0.7747 ; 

0.7366 

0.7435 

inital strength J 

SO 

0.4795 

0.4452 

0.4513 

mass density 

DMOD 

0.3134 

0.2938 

0.2968 

thickness 

THICK 

0.2298 

0.3889 

0.3722 

current temperature 

T 

0.0814 

0.0760 

0.0770 

ulimate temperature 

TU 

0.0795 

0.0740 

0.0749 

empirical constant 

QQ 

0.0571 

0.0534 

0.0540 

empirical constant 

W 

0.0517 

0.0860 

0.0491 

Poisson's ratio 

PMOD 

0.0064 

0.0040 

0.0041 

current number of cycles 

N 

0.0059 

0.0056 

0.0056 

reference temperature 

TO 

0.0036 

0.0033 

0.0034 

ultimate number of cycles 

NU 

0.0033 

0.0032 

0.0032 

reference number of cycles 

NO 

0.0026 

0.0049 

0.0024 

modulus of elasticity 

EMOD 

0.0021 

0.0019 

0.0041 
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Figure 1 .6. 1 below, is a graphical representation of the sensitivity factors results. 

What is important to note is that the four random variables with the highest sensitivity factors 
are the same for both the fully correlated and the partially correlated random fields. Therefore, 
the assumption of a fully correlated blade thickness is justified for future analysis. 


Sensitivity Factors for the Random Variables 


1.0 r 



THICK PMOD DLOAD N W TU TO 

EMOD DMOD NU NO SO T QQ 

Variables 


Figure 1.6.1 Sensitivity Factors for the Random Variables 
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1.7 DISTRIBUTION TYPE ANALYSIS 

The type of distribution for each random variable is based on experimental data, which 
is frequently lacking. There is some uncertainty in choosing the type of distribution to most 
accurately represent the data. Decisions are often made on the basis of tradition or experience. 
Choices are often controversial. However, given the mean and standard deviation, the 
following guidelines exist[16]. 

• Normal: Useful if COV is “small”, say (less than. 10%) 

- Tolerances 

- Modulus of elasticity 

- Poisson’s ratio 

- Various material properties 

• Lognormal: Good “default” distribution. Can be used for any variable. Plays 

important role in probabilistic design. 

- Cycles to failure in fatigue 

- Material strengths 

- Loading variables 

• Weibull: Very popular distribution, but probably overused. 

- Fatigue 

- Material strength 

- Time to failure in reliability analysis 

- Long-term distribution of stress ranges in fatigue 
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Probability 


To see the effect of distribution type on the output, several analyses were conducted. 
Initially, the distribution types for all random variables were assumed to be normal. Figure 
1.7.1 illustrates the effect distribution type had on the probability. 


Probability Density Plot for The Loading (DLOAD) 



Figure 1.7.1 Probability Density Plot for the Centrifugal Loading (DLOAD) 
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F.ffoct of Distribution Type op Probabilit y Results 

Table 1.7.1 summarizes the probability results for changes in distribution types for 
both fully and partially correlated random fields. Note that the reliability did not change 
significantly for a change in distribution type except for the change in distribution to Weibull 
for the centrifugal load. It has a reliability of 0.9999 or the so-called 4-9's of reliability. 


Table 1.7.1 Effect of Distribution on Probabilities of Failure 


Blade thickness 
uncertainty 

Random 

Variable 

Changed 

Distribution for 
Random 
Variable 
Changed 

Probability 
of Failure 
(AMV+) 

Reliability 
Based on 
AMV+ 

1 fully correlated 
random field 

none 

no change 

2.207E-04 

0.999779 

1 fully correlated 
random field 

Load 

lognormal 

3.590E-04 

0.999641 


Load 

Weibull 

1.386E-05 

0.999986 


Strength 

lognormal 

1.979E-04 

0.999802 


Strength 

Weibull 

5.678E-04 

0.999432 

1 fully correlated 
random field 

Density 

lognormal 

2.286E-04 

0.999771 


Density 

Weibull 

2.203E-04 

0.999780 


Thickness 

lognormal 

2.191E-04 

0.999781 

1 fully correlated 
random field 

Thickness 

Weibull 

2.384E-04 

0.999762 
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Effect of Distribution Type on Sensitivity Factors 

Initially, all analysis were performed with the 13 random variables and the thickness 
random field assumed to have normal or Gaussian distributions. Then the distribution type for 
one variable at a time was changed from the initial normal to a lognormal and then to a Weibull 
distribution. 

1) Centrifugal Loading (DLOAD) 

Figure 1 .7.2 illustrates the effect of distribution type on the centrifugal loading 
(DLOAD) random variable. Changing the distribution type to the lognormal distribution 
changed the sensitivity factors only nominally. The order of the highest to lowest of the four 
“significant” sensitivity factors did not change. 

However, a change in distribution type to Weibull resulted in a change in both the 
magnitude and the order of the “significant” sensitivity factors. The order of the first two 
factors, strength (SO) and centrifugal loading (DLOAD), changed. The magnitude of the 
“significant” sensitivity factors increased from 3 1 .5% to 43.6 % with the exception of the 
loading variable (DLOAD) which decreased from 0.7968 to 0.5396, a 32.3% decrease. This 
was the result of the shift in the Weibull distribution to the right and as was noted before, 
Weibull tends to be a better distribution for material strength. 


Effect of Distribution on Loading ( DLOAD ) 



EMOD DMOD NU NO SO T QQ 

Variables 


Figure 1.7.2 Effect of Distribution on Loading (DLOAD) 
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2) Strength (SO) 

Figure 1.7.3 illustrates the effect of distribution type on the strength (SO) random 
variable. Changing the distribution type to the lognormal changed the sensitivity factors only 
nominally. The order of the highest to the lowest of the four “significant” sensitivity factors did 
not change. However, a change in distribution type to Weibull resulted in a change in both the 
magnitude and the order of the “significant” sensitivity factors. The order of the first two 
factors, strength (SO) and centrifugal loading (DLOAD), changed. The magnitude of the 
“significant” sensitivity factors decreased from 20.2% to 24.3%. With exception of the initial 
strength variable (SO) which increased from 0.4566 to 0.73443, a substantial 60.8 % increase. 


Effect of Distribution on Strength (SO) 
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Figure 1.7.3 Effect of Distribution on Strength (SO) 
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3) Density (DMOD) 

Figure 1 .7.4 illustrates the effect of distribution type on the mass density (DMOD) 
random variable. Changing the distribution type to the lognormal distribution changed the 
sensitivity factors only nominally. The order of the highest to the lowest of the four 
“significant” sensitivity factors did not change. Changing the distribution type to Weibull 
resulted in a nominal change in the sensitivity factors. Most of the “significant” sensitivity 
factors increased slightly from 3.9% to 5.2%. Only the centrifugal loading decreased by 
2 . 8 %. 


Effect of Distribution on Density ( DMOD ) 
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Figure 1 .7.4 Effect of Distribution on Density (DMOD) 
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4) Thickness (THICK) 

The final variable tested for distribution type sensitivity was the blade thickness random 
field. Figure 1.7.5 illustrates the effect of distribution type on the blade thickness. Changing 
the distribution type to the lognormal distribution changed the sensitivity factors only 
nominally. The order of the highest to lowest of the four “significant” sensitivity factors did 
not change. A change in distribution type to Weibull, however resulted in changes in the 
“significant” sensitivity factors varied from a 4.3% decrease to a 3.4% increase. Only for the 
blade thickness (THICK) did the magnitude increase from 0.2128 to 0.2812, a 32.5% 
increase. 

The selection of distribution type does affect the sensitivity factors of the centrifugal 
loading and initial strength. The choice of Weibull for these two random variables increases 
the initial strength. This result was supported by the guidelines given in the beginning of this 
section. 


Effect of Distribution on Thickness ( THICK ) 
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Figure 1.7.5 Effect of Distribution on Thickness (THICK) 
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1.8 CONCLUSIONS 


In this paper, a probabilistic structural analysis of an SSME turbopump blade utilizing 
the program NESSUS is reported. First, probabilistic analysis yielded a statistical distribution 
for the critical von Mises stresses at the blade root due to uncertainty in the following variables: 
Blade thickness. Young’s modulus, Poisson’s ratio and mass density. Then NESSUS, with 
the newly implemented MSD model calibrated for Inconel 718 yielded a distribution of the 
material strength for the combination of high temperature (1000 °F) and high-cycle fatigue 
effects. With both stress and strength variabilities quantified, the structural reliability of the 
turbopump blade was determined using both the AMV+ and Monte Carlo Simulation methods. 

A comparison of the probability results obtained using the AMV + method to the well 
known and accepted Monte Carlo Simulation method yielded a percent difference of 14.7% for 
a partially correlated random field thickness and a percent difference of 7.7% for a fully 
correlated random field thickness. However, the corresponding 200,000 sample points used 
for the Monte Carlo Simulation method required six days of computational time versus three to 
four minutes for the AMV+ computation. Thus, the AMV+ method yielded accurate results 
with much greater efficiency. 

An important advantage of the reliability approach over the factor of safety approach is 
the ability to assess which random variables contribute most to the variability of the overall 
system. Sensitivity factors provided this insight which can be used for design optimization. 
Sensitivity factor results showed centrifugal loading and initial strength to be the most 
significant random variables for both a fully correlated and a partially correlated random field 
thickness. Mass density and thickness were third and fourth depending on correlation. The 
remaining ten random variables were one to two orders of magnitude less and therefore 
assumed insignificant in influencing the response. 

Reliability results and sensitivity factor results varied only slightly for a difference in 
random field thickness correlation. Hence, only a few of the random variables comprising the 
random field thickness are significant and account for the majority of the uncertainties in the 
random field. Thus, the assumption of a fully correlated blade thickness will not decrease the 
accuracy of results and is justified for future analysis. 

The effects of distribution type on both reliability and sensitivity factor results were 
studied. In both cases, a change in distribution from Gaussian to Lognormal resulted in only 
nominal changes. However, the change in distribution type from Gaussian to Weibull resulted 
in an increase in reliability for the centrifugal load variable from 0.99779 to 0.999986 or the 
so-called 4-9's of reliability. Also, changing distributions to Weibull resulted in a change in 
both the magnitude and the order of the ‘significant” sensitivity factors. Changing the loading 
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and strength variables’ distributions to Weibull swapped the order of the first two sensitivity 
factors. A Weibull distribution for strength (SO) increased the magnitude of its strength 
sensitivity factor from 0.4566 to 0.73443, a substantial 60.8% increase. 
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2.0 EDUCATION 


2.1 INTEGRATION OF RESEARCH AND EDUCATION OBJECTIVES 

It is well understood that when research is carried out in an academic setting that 
additional education benefits may occur beyond simply completing the research and delivering 
the results to the funding agency. These benefits do not automatically occur, but have a higher 
probability of occurring if the project Principal Investigator incorporates education objectives 
into the project and designs a carefully thought-out plan to achieve them. NASA and other 
agencies recognize this potentially effective strategy and have recently begun to request that 
proposals integrate within a single grant both research and education objectives. This two year 
NASA/UTSA Partnership Award is an example of such a project. 

A short-term education objective carried out during this Partnership Award was to 
expose undergraduate students to a NASA project in the hopes of encouraging them to obtain 
an advanced degree. To this end, undergraduate and graduate students were supported in 
faculty supervised research. Mr. Cody Godines was identified as a good and promising 
undergraduate student in mechanical engineering. He was recruited to join the first year 
Partnership Award project and was subsequently hired as an Undergraduate Research 
Assistant. As a result of working on the project he expressed an interest in continuing his 
education at the graduate level. He graduated in December 1998 with a BS in Mechanical 
Engineering and enrolled as a graduate engineering student during the second year of this 
grant. 


2.2 1998 AND 1999 SUMMER COURSES 

One long-term education objective of this project was to introduce undergraduate 
engineering students to probabilistic structural analysis and reliability via NESSUS. ME 4953 
Special Studies in Mechanical Engineering: Finite Element Applications in Solid Mechanics and 
Design, a FEM structural applications course with an introduction to NESSUS was conducted 
during the 1998 and 1999 summer semesters. These courses attracted over 65 mechanical 
engineering students. Approximately one third of the 65 students who completed these courses 
were minority students. 

A unique aspect of the course is that it was team taught. The instructor of record and 
lead laboratory instructor for the 1998 summer course was Ms. Callie C. Bast. Dr. Lola Boyce 
was the course lecturer and Mr. Mark T. Jurena was the assistant laboratory instructor. The 
laboratory instructors were budgeted from grant funds. Without the Partnership Award funds, 
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Ms. Bast would have had no assistance with the course and the quality of the course would 
have suffered. 

The instructor of record for the 1999 summer course was Dr. Randall D. Manteufel. 
Ms. Callie C. Bast and Mr. Cody Godines were the laboratory instructors. All three 
instructors for this summer course were supported by grant funds, and the course would not 
have been offered if not for this Partnership Award. In addition, this grant funded one 
engineer from SwRI, Mr. David Riha, to present guest lectures and assist with the computer 
laboratory portions of the summer courses. 

Both the 1998 and the 1999 summer courses consisted of a significant laboratory 
component utilizing student-centered learning activities. Student-centered learning includes 
laboratory-type assignments and projects on UNIX workstations wherein the students must 
directly and actively participate. This is in contrast to the lecture-type format wherein the 
learning is professor-centered and students participate via passive listening. This type of 
experiment in student-centered learning is receiving wide attention in Schools and Colleges of 
Engineering in the United States. Syllabi for both summer courses are provided in Appendix I. 

2.3 NESSUS STUDENT USER’S MANUAL 

A NESSUS Students User’s Manual was initiated during the first year of this 
Partnership Award and completed during the second year. It includes a brief overview of the 
program, explanations of the minimum number of NESSUS keywords necessary to work 
laboratory example problems, explanation of output files and a set of example 
problems/assignments drawn from structural analysis and reliability applications. This manual 
was used in both summer courses and is planned to be used in the follow-on courses scheduled 
for Spring 2000 and Summer 2000. This student manual is provided in Appendix II. 
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3.0 ACCOMPLISHMENTS 


The goal of this NASA Partnership Award was to advance innovative research and 
education objectives in theoretical and computational probabilistic structural analysis, 
reliability, and life prediction methods for improved aerospace and aircraft propulsion system 
components. This grant resulted in significant accomplishments in research and education, and 
the enhancement of UTSA’s engineering research environment. It allowed four UTSA 
Mechanical Engineering students; Mr. Mark T. Jurena , Mr. Cody Godines, Mr. Henock Perez 
and Mr. Ricardo Ramirez to work directly with the principal investigator, Ms. Callie C. Bast, 
providing them with a unique research experience that otherwise would not have been possible 
without this grant. 

3.1 ACCOMPLISHMENTS: RESEARCH 

Probabilistic structural analysis and reliability methods support the specific needs 
within the Aeronautics and Space Transportation Technology Enterprise by contributing to the 
development of a next-generation design tool that promises to improve the safety and reliability 
of both civil aviation and reusable space launch vehicles. For example, an enhanced NESSUS 
code is now available to assist in the design of the reusable Space Shuttle Main Engine (SSME) 
through the calculation of the reliability of turbopump blades. Thus, it is now possible to 
consider future engine designs based upon target reliabilities through a practical next-generation 
design tool. 

NESSUS Enhancement. Verification and Evaluation 

In order to provide NESSUS with the capability of determining the reliability of a 
structural component utilizing its material strength variability, the MSD model (supported by 
PROMISS, a probabilistic material strength degradation program developed by NASA Lewis 
Research Center with UTSA) was implemented in NESSUS in year 1 by SwRI, whom is 
responsible for maintaining the code. Thus, they assured that the new capabilities of NESSUS 
were fully integrated into the code and that they were compatible with previous capabilities. 
Initial verification of the enhanced NESSUS code (version 6.2), conducted in year 1, included 
two MSD model effects: high cycle mechanical fatigue and high temperature. The cumulative 
distribution functions (c.d.f.s) of lifetime strength were calculated using both PROMISS and 
NESSUS. Agreement of results was good in the lower tail but differed in the upper tail. 
Discrepancies were postulated to possibly be due to the difference in random number 
generators used in the two codes. It was determined that this phenomenon would be 
investigated in the proposed second-year continuation of this Partnership Award. Upon further 
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investigation in year 2, a sampling bias in the PROMISS code was found. The same starting 
seed was set for sampling each random variable. This bias was easily corrected in the 
PROMISS code. Verification of the enhanced NESSUS code was then completed. 

During year 2, this enhanced version of NESSUS was exercised and evaluated through 
a reliability analysis of an SSME turbopump blade (see section 1 of this report). This effort 
resulted in a thesis for one of the students supported on this grant 

Also in year 2, convergence criteria were improved for the most probable point based 
methods (MPP) in NESSUS, resulting in a new enhanced NESSUS code, version 6.3 (see 
Appendix III, SwRI Final Report). Southwest Research Institute also provided UTSA 
researchers and students with technical assistance in using the NESSUS code. Mr. David 
Riha, provided guest lectures and assisted with the computer laboratory for both the 1998 and 
1999 summer courses. The complete second year NESSUS Enhancement and Technology 
Support Final Report is provided in Appendix III of this report. 
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3.2 ACCOMPLISHMENTS: EDUCATION 

The introduction of undergraduate students to probabilistic finite element analysis 
methods was achieved through a UTSA summer course, ME 4953 Special Studies in 
Mechanical Engineering: Finite Element Applications in Solid Mechanics and Design. This 
FEM structural applications course with an introduction to NESSUS was conducted during the 
1998 and 1999 summers. Over 65 mechanical engineering students enrolled for this summer 
course. Students successfully completing it received three semester-hours of credit for a 
senior-level technical elective. A large number of the students who completed the course were 
minority students. 

A NESSUS Student User’s Manual was initiated in Year 1 and completed in Year 2. 
This manual will provide guidance in using NESSUS for future courses and help insure the 
continuation of probabilistic finite element analysis courses at UTSA. 

3.3 STUDENT ACHIEVEMENTS 

A graduate student, Mr. Mark Jurena, was recruited in the second-year continuation of 
this Partnership Award to complete his Master of Science in Mechanical Engineering degree. 
Supported by this Partnership Award, his thesis topic directly related to the research objectives 
of this grant. He presented his thesis “work in progress” at the AIAA/ASME/ASCE/AHS/ASC 
Structures, Structural Dynamics, and Materials (SDM) Conference last April in St. Louis, 
Missouri. Having concluded his research, he presented his research results in his thesis 
entitled “Structural Reliabilty Using Finite Element Analysis And A Probabilistic Material 
Strength Degradation Model”, and successfully completed his MSME degree in August, 1999. 

An undergraduate student, Mr. Cody Godines, was recruited in Year 1 of this 
Partnership Award as an Undergraduate Research Assistant. He earned his Bachelor of 
Science in Mechanical Engineering last December. While supported by this grant, Mr. 

Godines worked on a probabilistic finite element analysis of a beam and presented his 
NESSUS project at a NASA GRC-sponsored HBCU/HSU (Hispanic-Serving University) 
Research Conference in April, 1999. His conference paper is included in Appendix IV of this 
report. During this time, he became interested in going to graduate school. He enrolled as a 
graduate student in Spring 1999. Mr. Godines plans to complete his MSME degree in 2001 
under the guidance of Dr. Randall Manteufel and the support of a 1999 Partnership Award for 
Innovative and Unique Education and Research Projects. 
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The University of Texas at San Antonio 


Syllabus 

ME 4953 Special Studies in Mechanical Engineering: 
Finite Element Applications in Solid Mechanics and Design 
Summer, 1998 


INST RUCTORS: 


LECTURE/LARi 


Course Instructor of R ecord and Lead Laboratory Instructor 
Callie Bast, MSME; 

Office: 1.04.06 EB; 458-5588; Office Hours: TR, 2: 15 - 2:45 pm 
and 4:45 - 5:30 pm. 

Course Lecturer: 

Lola Boyce, Ph. D., P.E., 

Office: 3.04.42 EB; 458-5512; Office Hours: By appt. 

Assistant La boratory Instructor: 

Marie Jurena, BS 

Lecture: TR 2:45 pm to 3:40 pm 

Lab: TR 3:40 pm to 4:35 pm 


PREREQUISITES; 

TEXTBOOK; 

COMPUTER USAGE: 

COURSE OUTLINE; 


COURSE OBJECTIVE; 


Permission of Instructor. 

No required text. Course given from instructor’s notes, handouts, 
and recommended reading from various texts. 

ABAQUS, PATRAN, and NESSUS are available for assistance 
with projects, reports, and exams. 

A review of the necessary mathematical tools for the finite element 
method (FEM). Approximate numerical solutions to selected solid 
mechanics examples posed as FEM problems. Introduction to 
probabilistic methods in solid mechanics and mechanical design 
posed as FEM problems. 

To provide an understanding of analysis, design, and applications in 
solid mec hani cs and mechanical design using the finite element 
method (FEM) and appropriate FEM software tools. 


GRADING: A final grade will be assigned based on the 

following percentages: 


Midterm Exam 25% 

Final Exam 30% 

Homework/Labs 25% 

Project 20% 


EXAMINATION SCHEDULE: Midterm Exam: Thursday, June 25, 1998 

Final Exam: 3:15 - 5:45 pm; Tuesday, August 4, 1998 


ATTENDANCE: Mandatory. Attendance will be taken randomly. 

Note: Homework is assigned weekly and is good preparation for examinations. Questions regarding homework 
and other assignments are considered during office hours. No late assignments are accepted. Make-up exams are 
given only for extenuating circumstances discussed in advance and are generally more difficult than die regularly 
scheduled exam. 
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ME 4953 Special Studies: Finite Element Applications in Design 

Summer 1999 

Lecture: T-Th 2:45-4:35 EB 2.04.06 
Laboratory: T-Th 1:00-2:30 Session 001 
Laboratory: T-Th 5:00-6:30 Session 002 


Instructor: Randall D. Manteufel, Ph.D., P.E. 

Office: EB 3.04.54 

Phone: 458-5522 

Office Hours: T-Th 1 :30-2:30 


Lab Instructor: Callie Bast, MSME 

Office: EB 1.04.06 

Phone: 458-5588 

Lab Assistants: Mark Jurena, BSME 
Cody Godines, BSME 
Office: EB 1.04.06 

Phone: 458-5586 


Text: 


PWS-Kent, 1993. 


sign. C.E. Knight. 


Pro/MECHANICA Structure Tutorial, 2 ”* ed. Roger Toogood, SDC 
Publications, [http:// www.amazon.com /, ships in 24 hours] price $49.95. 

Description: A review of the mathematical foundations of the finite element method 

(FEM) for problems in continuum mechanics. Design case studies and finite 
element applications in mechanical design. Introduction to probabilistic 
FEM design. This class will emphasize the application of the FEM using 
commercial software. 


Course Goals: To give students the opportunity to learn and be able to: 

1) understand and explain the mathematical foundations of finite element 
analysis 

2) understand and explain numerical algorithms used in the finite element 
method 

3) understand and explain probabilistic methods used in mechanical design 
and reliability analyses 

4) perform mechanical analysis of components using finite elements 

5) use commercial software for finite element analysis in mechanical design 

Prerequisites: ME 3813 Solid Mechanics 

ME 3423 Applied Engineering Analysis 

Grading: Homework/Laboratory Assignments 50% 

Individual Design Project 50% 
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The University of Texas at San Antonio 


STUDENT MANUAL FOR THE PROBABILISTIC 
FINITE ELEMENT CODE, NESSUS 
(Numerical Evaluation of Stochastic Structures Under Stress) 



NASA GLENN RESEARCH CENTER 
PARTNERSHIP AWARD 
NAG # 3-2060 
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NESSUS 


Introduction 

Historically, NESSUS has provided probabilistic structural analysis for the response of 
engineered structures. It integrates both probabilistic analysis and finite element methods to 
produce a useful tool for both design and analysis of critical structural components. NESSUS 
allows for both random and deterministic descriptions of input and output quantities. For 
example, input parameters for material properties such as Young’s modulus, Poisson’s ratio, 
and mass density can be modeled as random variables wherein the user provides the 
appropriate statistics, namely mean, standard deviation and distribution type [2]. In addition, 
random input can include random fields, as well as random variables. Distributions such as 
Gaussian, Weibull, or lognormal may be selected. Other input parameters, such as loading 
type, geometry, boundary conditions and initial conditions can also be considered random. 
Then, probabilistic analysis yields the cumulative distribution functions (CDF’s) of the output 
random variables, e.g., stresses, strains, displacements, etc. 

Code Organization 

NESSUS is divided into ten parts. The preprocessor, NESSUS/PRE, provides for the 
representation of random fields. NESSUS/PRE computes a set of independent random 
variables from the random fields using an eigenvalue decomposition [3]. PFEM is the module 
for coordinating the operations between the finite element code and FPI. Next, a deterministic 
finite-element module, NESSUS/FEM, generates the deterministic responses of the finite- 
element structural model to a number of prescribed perturbations of the input variables. 
NESSUS/FEM contains a mixed iterative procedure to obtain the response of a structure. A 
Newton-Raphson iterative procedure based on a mixed variational principle can be used in 
NESSUS/FEM. Besides FEM, BEM, the boundary element module is used for structural 
analysis. These deterministic solutions provide the data from which a first or second-order 
Taylor series of the response is fit that defines the explicit functional relationship between the 
input random variables and the response of interest. The postprocessor, NESSUS/FPI (Fast 
Probability Integration), performs the probability calculations. The postprocessor also includes 
a choice of probabilistic analysis methods other than FPI, including for example, Adaptive 
Importance Sampling (AIS) and Monte Carlo simulation [4,5]. SIMFEM however, contains 
the driving module for the Monte Carlo as well as the Latin Hypercube sampling of the 
structure. There are also other packages for specialized analysis such as the RISK, SYSTEM 


N AS A/CR— 2001-211112 


57 



and SYSRSK modules. The RISK module computes the risk with respect to cost or a user 
defined criteria. The SYSTEM module is an alternate system reliability method for 
probabilistic fault tree analysis. The SYSRSK module computes the system risk by combining 
the cost of failure, in terms of replacement, inspection, repair or other criteria of individual 
components of the system with their probability of failure. The current version of NESSUS 
provides these ten component parts in an integrated package. These modules can be used in 
various combinations as the user desires. 


Input and Output Files 

NESSUS utilizes many different files. Filenames are based upon a logical naming 
convention. Depending upon input quantities within certain *PRINT or *PRINTOPTION 
statements, various output files will be created. The primary files are as follows: 


myfile.dat 

myfile.out 

myfile.pdb 

myfile.fem 

myfile.mov 

myfile.zal 

myfile.fpi 

myfile.fpo 

myfile.rvd 


This is the name of the input file. It contains information pertaining to finite 
element input, probabilistic input and Fast Probability Integration input. 

This file contains the majority of the output. It is the output most used for 
determining if the input data was correct. Explanations of errors are detailed 
within this file. 

This is the perturbation database. It is used by NESSUS to extract data for 
probabilistic analysis and is in a binary format. 

This is the deterministic input analyzed by the NESSUS/FEM module. 

This is the primary results from the probabilistic analysis. It contains 
information in a concise format for quick evaluation. 

A probabilistic results file suitable for input into a spreadsheet such as 
Microsoft Excel and subsequent graphical output. 

An intermediate output from NESSUS containing only the FPI input values. 
The following files are occasionally seen with various print statements. 

An intermediate output from NESSUS containing the probabilistic output. 

An intermediate output file form NESSUS. It is usually seen as the output 
from the NESSUS/PRE module. 


myfile.feo The deterministic output from the NESSUS/FEM module. 
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NESSUS Keywords 


♦PFEM The Probabilistic Finite Element Analysis Method 

This is the section that defines what random variables are input into NESSUS 
and what kind of analysis is performed on them. 

♦ZFDEFINE The beginning of the z- function or performance function definition. This 

defines what kind of response function is what adding variables that resist the 
analysis input. 

♦ZFUNCTION {Resistance model #} {# of coefficients} Selects the response function. 

{zcoef} Listing of real coefficients to be made available to the response function. 

C { Resistance model #} = 1 Structural design factor model. 

This predefined resistance model defines 
(strength - stress). 

= 2 Structural design factor model. 

This predefined resistance model defines 
(stress - strength). 

= 3 Material Strength Degradation (MSD) model. 
This predefined resistance model is explained in 
detail in a later section of this manual. 

♦MSDM Only required when using Resistance model #3. (See MSD model section of 
this manual) 

♦COMPUTATIONAL METHOD {method #} {number of random variables} 

C »> comments {method #} selects the type of analysis. 

C 1 corresponds to the Finite Element Analysis method 

C {number of random variables} : the number of computational random variables. 

{integer list of random variables} A list of the random variable numbers (12 3, etc.) 


♦END Signifies the end of the input for the ZFDEFINE section. 

♦RVDEFINE This is the start of random variable definitions. 


♦DEFINE {random variable #} Designation # of random variable. 

{random variable name} The name of the random variable. This is useful for the user to keep 

track of the random variables; give them meaningful names. 

{mean#} {standard deviation #} {distribution type} 

C {mean#} The mean of the distribution. 

C {standard deviation #} The standard deviation of the distribution. 


C 

C 

C 

C 

C 

C 

C 

c 

c 

c 

c 


{distribution type} The types of distributions in NESSUS: 

NORMAL -Normal or Gaussian Distribution 

WEIBULL 

EXTREMEVALUE 

LOGNORMAL 

CHISQUARE 

MAXENTROPY 

NESSUS 

FRECHET 

TWEIBULL Truncated Weibull Distribution 

TNORMAL Truncated Normal Distribution 
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{type} 

C 

C 

c 

c 

c 

c 

c 

c 

HARMONIC 

{databock} 


Defines the type of finite element data this random variable corresponds to. 
Includes the following types: 


ACCELERATION 

BEAMSECTION 

COEFFICIENTS 

COORDINATES 

DAMPING 

DISPLACEMENT 

DISTRIBUTED LOAD 

FORCES 


ORIENTATION 

PRESSURE 

PROPERTIES 

PSD 

SPRINGS 

TEMPERATURE 

UPERT 

VELOCITY 

YIELDFUNCTION 


Defines a block of additional data as required by an input. See expanded 
explanation in the *DEFINE {type} section. 


*PERTURB {perturbation #} Defines an individual perturbation 
{random variable #} {standard deviation shift} 


C {random variable #} The random variable that is to be perturbed. 

C {standard deviation shift} The number of standard deviations that an individual 

random variable will be perturbed or shifted. 


C MEAN VALUE PROBABILISTIC ANALYSIS 


*MVDEFINE This signals the start of the section where the keywords and random 

variables will be analyzed using the mean value first or second order 
methods. States type of data to be analyzed and which nodes or 
components are to be analyzed. 


*PERT {# of perturbations} Selects the total # of perturbations to be used in the 

probabilistic analysis. 

{list of perturbations} an integer list of perturbations to be perturbed 


* RAN VAR {# of random variables} Selects the total number of random variables to be analyzed, 
{random variable # list} an integer list of random variable numbers to be perturbed 


*DATATYPE {type of data #} 

C {type of data #} 

C 

C 


Specifies the type of data on which to perform the 
probabilistic analysis. 

= 0 specifies Incremental 
= 1 specifies Eigenvalue 
= 2 specifies Harmonic/spectral 
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*RESPTYPE {response type #} The type of response variable to extract from the 
perturbation database. 

C {response type #} 


C 1= TOTAL DISPLACEMENT 

C 2 = TOTAL STRAIN 

C 3 = TOTAL STRESS 

C 1 1 =PLASTIC STRAIN 

C 12 = BACKSTRESS 

C 13 = CREEP STRAIN 

C 14 = THERMAL STRAIN 

C 17 = GENERALIZED STRAIN 

C 18 = GENERALIZED STRESS 

C 21= MATERIAL STATE 

C VARIABLES 

C 25 = VELOCITIES 

C 26 = ACCELERATIONS 

C 30 = THE EIGENVALUE FOR 

C THE MODE 

C 31= MODAL DISPLACEMENT 

C (EIGENVECTOR) 

C 32 = MODAL STRAIN 

C 33 = MODAL STRESS 

C 35 = THE FREQUENCY IN 

C RADIAN PER TIME 


36 = THE FREQUENCY IN CYCLES PER TIME 

51 = REAL DISPLACEMENT 

52 = REAL STRAIN 

53 = REAL STRESS 

61 = IMAGINARY DISPLACEMENT 

62 = IMAGINARY STRAIN 

63 = IMAGINARY STRESS 

71 = THE AMPLITUDE OF THE DISPLACEMENT 

72 = THE AMPLITUDE OF THE STRAIN 

73 = THE AMPLITUDE OF THE STRESS 

8 1 = THE PHASE OF THE DISPLACEMENT 

82 = THE PHASE OF THE STRAIN 

83 = TIIE PHASE OF THE STRESS 

91 = MEAN SQUARE DISPLACEMENT 

92 = MEAN SQUARE STRAIN 

93 = MEAN SQUARE STRESS 
96 = STRESS VELOCITY 


♦CONDITION {beginning condition #} {ending condition #} 


Selects the beginning and ending 
condition numbers for the mean 


value analysis. 

C Select condition type with previous DATATYPE keyword. The condition refers to 

C here refers to either the increment number, the mode number or the 

C harmonic/spectral case number respectively. 


♦NODE {beginning node #} {ending node #} Selects the beginning and ending node 

numbers for the mean value analysis. 

If there is only one node, then only one 
number is needed. 


♦COMPONENT {beginning component #} {ending component #} 


Selects the beginning and 
ending component numbers 
for the mean value analysis. 


C 

C 

C 

C 

C 

C 

C 


Component numbers: 

1 - translation in the x direction 

2 = translation in the y direction 

3 = translation in the z direction 

4 - rotation with respect to the positive x direction 

5 = rotation with respect to the positive y direction. 

6 = rotation with respect to the positive z direction 
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* PRINT {print control} Controls the amount of printed output that will be generated 

during a PFEM analysis. 

All finite element and fast probability output sent to 
the .feo and .fpo files, respectively. 

All finite element and fast probability output sent to 
the .feo and .fpo files, respectively 
All output is sent to the .out file. 

*END Signifies the end of the MVDEFINE section. 

C ADVANCED MEAN VALUE PROBABILISTIC ANALYSIS 
♦AMVDEFINE Signifies the start of the advance mean value analysis definition. 

♦ITERATION Defines the convergence criteria for the AMV+ z-level and p-level iteration 

procedures. 

{maximum # of iterations} {error tolerance} 

{ maximum # of iterations } : The maximum number number of iterations 

the procedure is allowed to perform before 
exiting 

{error tolerance} : The relative error between successive 

iterations. 

Using the Zlevel Algorithm, convergence is defined to be 

= IP,- Pj/Pi-i = e when using the relative change in beta and 

= I Z,-Z, 1 |/Z,,= 0.2 when using the relative change in the Z 

and the angle between the MPP’s from the last two iterations, 0, is less 
than the allowable maximum of 30°. 

Theta is defined as cos 0= a ( • a^, , where a, are the direction cosines to 
the MPP from integration i . 

Using the Plevel Algorithm, convergence is defined as 
= | Zj - Z M |/ Z M = £ and 

the angle between the MPP’s from the last two iterations, 0, is less than the 
allowable maximum of 30°. 

Theta is defined as cos 0= a, • a M , where a, are the direction cosines to 
the MPP from integration i 

♦CONDITION {beginning condition #} {ending condition #} Selects the beginning and ending 

condition numbers for the 
Advanced Mean Value analysis. 

♦NODE {beginning node #} {ending node #} Selects the beginning and ending node 

numbers for the Advanced Mean Value 
analysis. 

♦COMPONENT {beginning component #} {ending component #} Selects the beginning and 

ending component numbers 
for the Advanced Mean Value 
analysis. 


C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


C {print control}: = SHORT 

C 

C = MEDIUM 

C = LONG 


♦END Signifies the end of the AMVDEFINE section. 
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♦END 


End of PFEM input. 

♦FEM Signals the start of the deterministic Finite Element section. 

♦CONS {constitutive model #} Selects the constitutive material model. 

C {constitutive model #} = 0 for linear elastic model. 

C = 1 selects the simplified plasticity 

C =2 selects classical von Mises J2-flow plasticity 

C =3 selects Walker's creep-plasticity model. 


♦ELEMENTS {# of elements} Allocates memory for a maximum number of elements, 

{element type #} The element type used for analysis. 

The following element type #’s are valid: 
bilinear, isoparametric plane stress element. 

trilinear isoparametric, solid element, 
bilinear, isoparametric, axisymmetric element, 
bilinear isoparametric plane strain element, 
bilinear, isoparametric, variable-thickness shell element, 
linear, isoparametric, tapered beam element, 
bilinear, isoparametric, plane stress element, 
bilinear, isoparametric, plane strain element, 
bilinear, isoparametric, axissymmetric element, 
trilinear, isoparametric, solid element 

♦NODES {# of nodes} Allocates memory for the number of nodes on the mesh. 

♦BOUNDARY {# of boundary conditions} Allocates memory for the boundary conditions. 


C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 


{element type #} 

3 = 

7 = 

10 = 

11 = 

75 = 

98 = 

151 = 

152 = 

153 = 

154 = 


♦FORCES {# of forces} Allocates memory for the applied forces. 

♦PRINT {# quantities printed} This controls the amount of data printed to the output file. 

C {# quantities printed} The maximum number of quantities printed to the output 

C file. The default is 10. 


♦MONITOR {# of monitored values} Monitors a small number of values when the code is 

run interactively. This shows what is happening to 
a certain quantity. 


♦END The end of the memory allocation. 

♦ITER {begining* perturbation #} {ending perturbation #} 

{maximum# of perturbations} {relative error} {absolute error} {displacement error} {energy error} 
C This sets the convergence tolerances for the iterative solution algorithms. 
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♦MONITOR To monitors certain values at the end of an iteration. 

{Monitored quantity} {Location} {Component#} This selects a small number of values that 

will be reported at the end of each iteration. 
Used primarily when running the code 
interactively. 


C 

C 

C 

C 

C 

C 

C 

C 


{Monitored quantity} The string denoting the quantity to be monitored. The following 
quantities are valid: 


CREEPSTRAIN 

EQUIVALENTPLASTICSTRAIN 

FORCE 

INCREMENT ALDISPLACEMENT 
PLASTICSTRAIN 
REACTION/RESIDUAL 


STRAIN 

STRESS 

TEMPERATURE 

THERMALSTRAIN 

TOTALDISPLACEMENT 


C {Location} The component or location to be monitored, i.e., NODE # 

C {Component#} The string specifying a particular component or location. The valid 

componets are: 

COMPONENT# 

NODE# 


♦COORDINATES Used to specify nodal point coordinates for the finite element mesh used in 

the analysis. 

{node#} {1 st-coordinate#} {2nd-coordinate #} {3rd-cooidinate #} . . . 


C 

C 

c 

c 

c 

c 


c 


For the element type 3,11 151 and 152, the two coordinates required are: 
{coordinate #}= 1 is the x global coordinate. 

2 is the y global coordinate. 

For the elements of type 10 or 1 53, the two coordinates are: 

{coordinate #}= 1 is the (z) axial coordinate. 

2 is the (r) radial coordinate. 


For the elements of type 
{coordinate #}= 1 is 

2 is 

3 is 


7 or 154, the three coordinates are: 
the x global coordinate, 
he y global coordinate, 
the z global coordinate. 


For the three-dimensional shell element type 75, the four coordinates are: 
{coordinate #}= 1 is the x global coordinate. 

2 is the y global coordinate. 

3 is the z global coordinate. 

4 is the thickness of the shell. 

This defaults to the last nonzero thickness entered. 

For the three-dimensional meshes using element type 98, the six coordinates are 
{coordinate #}= 1 is the x global coordinate. 

2 is the y global coordinate. 

3 is the z global coordinate 

4 is the x-component of the beam normal 

5 is the y-component of the beam normal 

6 is the z-component of the beam normal. 
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♦ELEMENTS {element type #} 

{element #} {node # a} {node # b} {node # c} {node # d} 

C These lines express element connectivity. 

C {element type #} The element type used for analysis. 

(See *ELEMENTS in *FEM section for valid #’s.) 

The total number of elements n must be less than or equal to the maximum specified in parameter 
data using the keyword option *ELEMENTS. The exact number m of nodes used to define the 
connectivity list for a given element will depend on the element type used for the analysis. It is important 
that the element connectivity list be entered in the correct sequence. For a two-dimensional, four-noded 
element, this means that the nodes must be entered in counterclockwise order, as shown in Figure 1. Th< 
node numbering convention for a three- dimensional, eight-noded element is also shown. The correct 
node order is as indicated by the alphabetic sequence A, B, C ... Z. If the node sequence for an elemeni 
is found to be incorrect, the code will attempt to fix it by reordering the nodes. A warning message will 1 
printed to the output file, listing the modified node order. When this happens, the user should always 
check the modified node sequence to see whether it still agrees with the desired mesh topology. 
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c 



Plane stress, 
plane strain, or 
axisymmetric 



Node numbering for continuum-type elements. 


•BOUNDARY To specify displacement boundary conditions. 

{node #} { component #} {prescribed displacement} 

{node #} The node number. 

{component #} The components of the boundary conditions. 

1 — translation in the x direction 

2 = translation in they direction 

3 = translation in the z direction 

4 = rotation with respect to the positive x direction 

5 = rotation with respect to the positive y direction. 

6 = rotation with respect to the positive z direction 
{prescribed displacement} The displacement prescribed for that particular direction at 

that particular node. 


•PROPERTIES {element type #} This defines the material properties for a particular set of nodes, 
{beginning node #} {ending node #} {prop #1 } {prop #2} {prop #3} {prop #4} {prop #5} {prop # 6} 


C 

C 

c 

c 

c 

c 

c 

c 

c 


{beginning node #} 
{ending node #} 
{prop #1} 

{prop #2} 

{prop #3} 

{prop #4} 

{prop #5} 

{prop #6} 


Specifies the first node number in a series. 

Specifies die last node number in a series. 

Specifes the element thickness for a plane stress analysis. 

A dummy paramenter otherwise. 

Specifies the elastic modulus in units of force/unit length.. 
Specifies Posssion’s ratio a dimensionless quantity. 

Specifies the thermal expansion coefficient in units of 
Length/temperature. Also called alpha. 

Specifies the mass density in unit of mass/length A 3 
Specifies the shear modulus for a plane stress element and 
Axisymmetric problems and three dimensional solid analysis, in 
units of force/1 ength A 2. Unused otherwise. 
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♦FORCES To specify applied forces. 

{node #} {component #} {force} Specifies the location, component and force applied to a 

node. 


C {node #} Specifies the node number 

C {component #} Specifies the direction of the force (see *BOUND ARY) 

C {force} Specifies the magnitude of the load applied . 


♦PRINT 


{quantity type} {beginning #} {ending #} {location parameter} {beginning #} {ending #} 

Controls the amount of data printed to the output file. 

The beginning and ending #’s are optional, otherwise all 
components will be printed. 

{quantity type} Specifies particular quantity to be printed. 

The possible quantity types are: 


CREEPSTRAIN 

EQUIVALENTPLASTICSTRAIN 

FORCE 

INCREMENT ALDISPLACEMENT 

PLASTICSTRAIN 

REACTION/RESIDUAL 

STRAIN 

STRESS 

TEMPERATURE 
THERMALSTRA1N 
TOT ALDISPLACEMENT 


{location parameter} 

The possible location parameters are: 

ELEMENTS 

1NTEGRATIONPOINTS 
LAYERS 
NODES 

♦END Specifies the end of all model data input. 

♦FPI Specifies the start of the Fast Probability analysis input setion. 

♦RVNUM {random variable #} The number of random variables to be analyzed by Fast 

Probability Integration. 
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*GFUNCTION {function type #} Defines the g-function approximation. 


C {function type #} 

C 
C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

♦DATASETNM {# of data sets} 


The type of function desired by the user. 

= 1 for a linear g-function approximation. 

N+l datasets are required for linear, where N is the number 
of random variables or pertubations. 

= 2 for a quadratic g-function approximation. 

2*N+1 data sets are required for the quadratic 
approximation. 

= 5 for a g-function of the form (strength-stress). Used for 
reliability analysis with the strength-stress form, i.e. R-S. 
Only the probability of failure, i.e. P[g(x)<0], will be 
determined. 

= 6,7,8 for a user defined response function. 
g-fUnction must be programmed in routine RESPON or 
USERES, or defined in the input deck provided that the 
g-function has a continuous first derivative. 

Designates the number of data sets for a particular problem. 
(See *GFUNCTION for the number required.) 


*METHOD {method type #} Defines the solution type used in FPI analysis. 


C 

C 

C 

C 

C 

C 

C 

C 

C 


C 

C 

C 

C 

C 

C 

C 

C 

C 

C 


C 

C 

c 

c 

c 

c 


{method type #} 


= 0 First order reliability method (FORM). 

= 1 Advanced first order reliability method (FPI). 

Computes both first order reliability method and 
advanced first order reliability method solutions. 
Parameter = 1 is recommended. 

= 2 Fast convolution in the regular-space (CONVX). 
FORM solution computed first. 

= 3 Fast convolution in the standard normal-space 
(CONVU). 

FORM solution computed first. 

= 4 Second order reliability method (SORM). 

FORM solution computed first. 

= 5 Importance sampling method with radius reduction 
factor (ISAMF). *MONTE keyword and data are 
expected in the Model section of FPI. 

Monte Carlo can only be used for user-defined 
response levels (ZLEVELS), *ANALTYP=1 . 

= 6 Conventional Monte Carlo method (MONTE). 

* MONTE keyword and data are expected in Model 
section. Monte Carlo can only be used for 

* ANALTYP = 0 or 1 . 

= 7 Same as parameter = 5 except that the radius is user- 
supplied (ISAMR).Limited to * ANALTYP = 1 . 

= 8 Adaptive importance sampling (linear 

surface)(AISl). "TTOL keyword and data are 
expected in Model section. 

Limited to * ANALTYP = 1 . 
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c 

c 

c 

c 


= 9 Adaptive importance sampling (quadratic surface) 
(AIS2). *ITOL keyword and data are expected in 
Model section. 

Limited to * ANALTYP - 1 . 


C 


= 1 0 Mean value method (MV). 


C 


= 1 1 Advanced mean value method (AMV). 


C 

C 

C 


= 1 2 Advanced mean value method plus (AMV +). 

*ITER keyword and data are expected in Model 
section. 


♦PRINTOPT {printout type #} 


C 

C 

c 


{printout type #} Defines the amount of data to be printed out from FPI. 
= 0 for the short print out. 

= 0 for long print out. 


♦ANALYTYPE {analysis type #} 

C {analysis type #} Specifies the analysis type. 


C 

c 

c 

c 

c 

c 

c 

c 

c 

c 


The types used are: 

= 0 FPI defined probability levels are used. These are 

automatic and are typically from —4.5 to 4.5 standard 
deviations from the mean of the response. 

= 1 User defined response levels(Z levels). The *ZLEVELS 
keyword and the Zlevels are required in the FPI model 
section. This is the most efficient method. 

= 2 User defined probability levies. The PLEVELS keyword 
and the probability levels are required in the FPI model 
section. 


♦END Signifies the end of the parameter data. What follows this first *END is 

the FPI model data section. 


♦ZLEVELS {# of levels} Defines the number of response levels; maximum =20 

♦PLEVELS {# of levels} Defines the number of probability levels. 


♦MONTE 

{# of samples} {seed#} {beta factor} 

C Specifies the data necessary for running Monte Carlo Simulations. This 

C line is required when Method=5, 6 or 7 in the FPI parameter section. 


C 

C 


C 

C 

C 


{# of samples} The number of random samples for each response level using * METHODS 
5 and 7. 

It is the total number of random sampels for METHOD 6. 

{seed #} The random number generator seed, however, it must be an integer. 

{beta factor} = 5 This defines the reduction factor for the sampling radius based on 

the minimum distance. 

= 6 This value is zero . 

= 7 This value is a user-supplied radius. 


*itol 

♦ITER 

♦CONFINTVL 

♦END End of input file. 
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Material Strength Degradation (MSP) Model 


The MSD model in the form of a postulated randomized multi-factor equation provides for 
quantification of uncertainty in the lifetime strength of components subjected to a number of 
diverse random effects. The mathematical expression for this MSD model is: 



A. tt -A. 

iU l 

A iU~ A iO 


a i 


( 1 ) 


where, A,, and Aj 0 are the current, ultimate and reference values, respectively, of a particular 
effect; aj is the value of the calibrated empirical material constant for the i* effect terms of the 

variables in the model; and S and S 0 are the current and reference values of material strength. 

Each term has the property that if the current value equals the ultimate value, the lifetime strength 
will be zero. Also, if the current value equals the reference value, the term equals one and 
strength is not affected by that value. The product form of equation ( 1 ) assumes independence 
between individual effects. This equation may be viewed as a solution to a separable partial 
differential equation in the variables with the further limitation or approximation that a single set of 
separation constants, a,, can adequately model the material properties. 

The model is currently set up for five fixed effects and 12 user defined effects that 
typically reduce lifetime strength. These effects are listed below: 

• High temperature 

• High-cycle mechanical fatigue 

• Low-cycle mechanical fatigue 

• Creep 

• Thermal fatigue 

• Up to 12 user defined (General) effects 


When expanded for all effects, the equation is as follows: 
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Temperature HCF LCF Creep Thermal Fatigue General 


where T is temperature, N is cycles, t is time, and A is a user-defined effect. The lower case 
exponents are the empirical material constants for each effect. The subscript, u, refers to an 
ultimate value, the subscript, o, is the reference value, and the non-subscripted term is the current 
value for the effect. 
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To increase model sensitivity for certain effects, a logarithmic base ten transformation is 
introduced resulting in the following form. 


S = S„ 


LOGjAjy ) - LOGjA ^ ) 
LOG{A il} ) - LOG(A i0 ) 


(3) 


The nature of the material data will dictate the use of the log transformation. The effects to be 
used and the model type (log transformation) are defined in the Z-function definition section of the 
PFEM input file. Any combination of effects can be selected and any of the terms can be 
considered random if desired. The form of the g-function used by NESSUS is 


zg = S 0 


n 



(4) 


where 

S G is the reference value of material strength 

Au is the ultimate value of the particular effect 

Ajo is the reference value of the particular effect 

A, is the current value of the particular effect 

a i is the empirical material constant for the particular effect 

a is the structural response as calculated by NESSUS/FEM, i.e. stress. 
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As an example, the equation is expanded for the two effects of high temperature and high- 
cycle mechanical fatigue as shown in equation (5) below: 


s 

1 7 “ 

1 

£ 

1 

<7 

’ LOG ( N v ) - LOG (N) ' 

S 0 ~ 

1 

1 

1 


LOG (N v ) - LOG (N a ) 


The MSD model was calibrated for INCONEL 7 1 8 by appropriate curve-fitted least squares linear 
regression of experimental data. Linear regression of the data for each effect resulted in estimates 
for the empirical material constants, as given by the slope of the linear fit. These estimates, 
together with ultimate and reference values, were used to calibrate the model specifically for 
Inconel 7 1 8. Lifetime material strength results, in the form of cumulative distribution functions 
(CDF’s), illustrate the sensitivity of lifetime strength to any one, or a combination of, the effects. 


The formulation defined in NESSUS is 


\T U -T 1 


’ LOG (N v ) - LOG (N) ' 

» 

|K? 

1 

!b? 


LOG ( N v ) - LOG ( N n ) 


By defining S to be 1.0 and not including a computational model to compute stress (o = 0.0), 
then NESSUS is solving the problem for S/So ■ The random variables are defined in Table 1 . 

Table 1. Random Variable Definitions 


Effect 

Symbol 

Distribution 

Mean 

Standard 

Deviation 

High Cycle 

Nu 

Normal 

1. Ox 10 |U cycles 

1 .Ox 1 0 V cycles 

Mechanical 

N 

Normal 

2.5x10 s cycles 

2.5x1 0 4 cycles 

Fatigue 

No 

Normal 

0.25 cycles 

0.025 cycles 

(at 75 °F) 

W 

Normal 

0.3785 

0.0114 

High 

Tu 

Normal 

2369.0° 

236.90° 

Temperature 

T 

Normal 

1000.0° 

100.0° 

(at 1000 °F) 

To 

Normal 

75.0° 

7.5° 


Q 

Normal 

0.2422 

■EEsUH 
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The relevant NESSUS input for this problem is shown in Figure 1. Note that to consider a term 
in an effect random, the alphanumeric name must match the name in the RVDEFINE input. An 
example of this is the highlighted NU variable shown in both the ZFDEFINE section and the 
RVDEFINE section in Figure 1. 


, o 

LINEAR 
LOG NU 


* ZFDEFINE 
♦EXPLICITMODEL 
1,2, 3,4,5, 6,7,8 
♦ZFUNCT 3 

♦MSDM 1.0 (S 0 

TEMP 
HCF 
END 

♦END 

♦RVDEFINE 
♦DEFINE 1 

NU 

1.0E10 1.0E9 
* DEFINE 2 
N 

2.5E5 2.5E4 

♦DEFINE 3 
NO 

0.25 0 

♦DEFINE 

w 

0.3785 
♦DEFINE 
TU 

2369.0 
♦DEFINE 
T 

1000.0 
♦DEFINE 
TO 

75.0 7 

♦DEFINE 8 

QQ 

0.2422 


025 

4 


236. 

6 

100 . 

7 

50 


TU 

N 


T 

NO 


TO 


NORMAL 


NORMAL 


NORMAL 


0114 NORMAL 


NORMAL 


NORMAL 


NORMAL 


0.0088 NORMAL 


QQ 


Figure 1. NESSUS Input for Example Problem 
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♦DEFINE (type) 


ACCELERATION 
nodel daccl 1 dacc21 daccml 
node2 daccl 2 dacc22 daccm2 
noden daccl n dacc2n ... daccmn 

where the values dace denote relative perturbations to the m components of the initial accelerations 
at each node. This means that, when die associated random variable is perturbed by one full 
standard deviation, a total change of dstdev times dacc2n should be observed in the second 
component of the initial acceleration for node n. The same general rule also holds for other 
components 


COORDINATES 

nodel dcoorl 1 dcoor21 ... dcoorml 
node2 dcoorl 2 dcoor22 ... dcoorm2 
noden dcoorl n dcoor2n ... dcoormn 

where the values dcoor denote relative perturbations to the m components of the nodal coordinates 
at each mesh point. The input format for DataBlock is identical to the one used to specify 
♦COORDINATES as discussed elsewhere in this section. 


The input format used for defining a perturbation variable associated with the Rayleigh damping 
coefficients is 
DAMPING 1 
dalpha dbeta 

where jdamp is set to 1, for Rayleigh damping, and dalpha and dbeta are the relative perturbations 
to the two Rayleigh damping parameters 

DAMPING 2 
imodel jmodel dratiol 
imode2 jmode2 dratio2 
imoden jmoden dration 

where jdamp is set to 2, for modal viscous damping, and the values of dratio denote the relative 
perturbations to the modal damping ratios. 

DAMPING 3 
imodel jmodel dratiol 
imode2 jmode2 dratio2 
imoden jmoden dration 

where jdamp is set to 3, for structural damping, and the values of dratio are again used to denote 
the relative perturbations to the modal damping ratios. In all three cases, the formats are identical 
to the ones used to prescribe a given damping model using the *DAMPING keyword option. 
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DISPLACEMENT 
nodel ddisl 1 ddis21 ... ddisml 
node2 ddisl 2 ddis22 ... ddism2 
noden ddisl n ddis2n ... ddismn 

where the values ddis are relative perturbations to the m components of the initial displacements at 
each node. 


DISTRIBUTEDLOAD jetyp 
ieleml jeleml indexl ddistl 1 ddist21 ... ddistml 
ielem2 jelem2 index2 ddistl 2 ddist22 ... ddistm2 
ielemn jelemn index3 ddistl n ddist2n ... ddistmn 

where the values ddist are relative perturbations to the m distributed loads for each element 
between ielem and jelem. Notice the flag jetyp, which is used to indicate the element type used for 
the analysis. 

FORCES 

nodel ndof 1 dforcel 
node2 ndof2 dforce2 
noden ndofh dforcen 

where dforce are the relative perturbations to the force (or moment) acting at degree- of-free dom 
ndof of node node. It is assumed that a set of unperturbed loads has been defined at these nodes 
and degrees-of-freedom (even though these unperturbed loads may have zero magnitude). 

The three different perturbation variables associated with harmonic excitation parameters are 
discussed next. 

HARMONIC jharm 

dfreq 

FORCES 

nodel ndof 1 damplil dphasel 
node2 ndof2 dampQ dphase2 
noden ndofh damplin dphasen 

where dampli and dphase are relative perturbations to the amplitude and phase of the harmonic 
loads acting at degree-of-freedom ndof of node node, and dfreq - is the perturbation to the 
excitation frequency. It is assumed that unperturbed harmonic force excitations have been 
previously specified at these nodes and degrees-of-freedom for harmonic case jharm. 

To define a perturbation variable associated with harmonic base acceleration parameters, the input 
format is 

HARMONIC jharm 
dfreq 

ACCELERATION 
nodel ndof 1 damplil dphasel 
node2 ndo£2 dampli2 dphase2 
noden ndofh damplin dphasen 

where dampli and dphase are the relative perturbations to the amplitude and phase of the harmonic 
base accelerations for degree-of-freedom ndof of node node, and dfreq is the perturbation to the 
excitation frequency. It is assumed that unperturbed base accelerations have been specified 
previously at these nodes and degrees-of-freedom for harmonic case jharm. 
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Finally, the input format used for defining a perturbation variable associated with the harmonic 
pressure loading parameters is 
HARMONIC jharm 
dfreq 

PRESSURE 

inodel jnodel damplil dphasel 
inode2 jnode2 damplQ dphase2 

inoden jnoden damplin dphasen 

where dampU and dphase are the relative perturbations to the amplitude and phase of the harmonic 
pressure loading defined at nodes inode through jnode, and dfreq is the perturbation to the 
excitation frequency. 


ORIENTATION 

inodel jnodel dalphal dbetal dgammal 
inode2 jnode2 dalpha2 dbeta2 dgamma2 

inoden jnoden dalphan dbetan dgamman 

where dalpha, dbeta, and dgamma denote relative perturbations to the three direction angles 
defining the material orientation in three dimensions. For two-dimensional analysis only the first 
angle needs to be specified. 


PRESSURE 
inodel jnodel dpressl 
inode2 jnode2 dpress2 

inoden jnoden dpressn 

where the values of dpress are relative perturbations to the nodal pressures acting at nodes inode 
through jnode. 


PROPERTIES jetyp 

inodel jnodel dpropl 1 dprop21 ... dpropml 
inode2 jnode2 dpropl 2 dprop22 ... dpropm2 

inoden jnoden dpropln dprop2n ... dpropmn 

where the values of dprop are used to denote relative perturbations to the m material properties 
assigned to nodes inode through jnode. Notice the flag jetyp, which is used to indicate the element 
type used for the analysis. 
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The three different types of perturbation variables associated with power spectrum excitations are 
discussed next. The input format used for defining a perturbation variable involving power spectra 
for nodal forces is: 

PSD jpsd 
df req 1 dpsdl 
dfreq2 dpsd2 
dffeqm dpsdm 
FORCES 

nodel ndof 1 damplil 
node2 ndof2 dampV 
nodenndofii damplin 

where dffeq and dpsd define changes to the profile of the power spectrum excitation curve, and 
dampli is the relative change in amplitude of the power spectrum excitation for degree-of- freedom 
ndof at node node. It is assumed that unperturbed power spectrum forces have been previously 
specified at these nodes and degrees-of-freedom for the jpsd power spectrum case. 

The input format used for defining a perturbation variable associated with the power spectrum 

excitation parameters for base acceleration is 

PSD jpsd 

df req 1 dpsdl 

dfreq2 dpsd2 

dffeqm dpsdm 
ACCELERATION 
nodel ndof 1 damplil 
node2 ndof2 dampdi2 

noden ndofii damplin 

where dffeq and dpsd define changes to the profile of the power spectrum excitation curve, and 
dampli is the relative change in amplitude of the power spectrum excitation for degree-of-ffeedom 
ndof at node node. Again, it is assumed that unperturbed power spectrum accelerations have been 
specified at these nodes and degrees-of-freedom for the jpsd power spectrum case. 

To define a perturbation variable involving the power spectrum excitation parameters for nodal 
pressures, the input format is 

PSD jpsd 
df req 1 dpsdl 
dfreq2 dpsd2 

dffeqm dpsdm 
PRESSURE 
inodel jnodel damplil 
inode2 jnode2 dampi2 

inoden jnoden damplin 

where dffeq and dpsd define changes to the profile of the power spectrum excitation curve, and 
dampli denotes the relative change in the amplitude of the random pressure loading acting at nodes 
inode through jnode. 
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In all cases, the number of lines defining the profile of the power spectrum excitation curve is 
limited by the parameter npsdp, and the number of data lines following FORCES, 
ACCELERATION, or PRESSURE is limited by the parameter npsep, both defined by the *PSD 
option in the parameter data block for the problem.. 


Perturbation variables involving ground spring stiffnesses may be defined using the following 

input format 

■"DEFINE jpvar 

dmean dstdev 

SPRINGS 

nodel ndofl dstiff 1 

node2 ndof2 dstiff2 

noden ndof n dstiffn 

where dstiff is used to denote the relative perturbation to the ground spring stiffness at 
degree-of-freedom ndof of node node. It is assumed that a set of unperturbed ground springs has 
been defined at these nodes and degrees-offreedom. 

The input format used for defining a perturbation variable associated with the nodal temperatures 
is 


TEMPERATURE 

inodel jnodel dtempl 1 dtemp21 . . . dtempml 
inode2 jnode2 dtempl 2 dtemp22 ... dtempm2 

inoden jnoden dtempln dtemp2n dtempmn 

where the values dtemp denote relative perturbations to the m temperatures at each mesh point. 

The input for most elements requires only one temperature value at each mesh point. . However, 
elements with multiple layers, such as the bilinear shell element (type 75), will need to have 
temperatures specified for each one of the m layers at a node. 

The NESSUS finite element code also provides a facility for defining general loading 
time-histories that are a function of one or more random variables. These random variables are not 
necessarily, associated with any particular loading type, and may even represent quantities that are 
not explicitly defined in the finite element code, such as engine power level, flow rate, and inlet 
temperatures in a duct, etc. The functional dependence of the loading on these random variables 
may be coded as a general analytical (or algorithmic) expression in subroutine UPERT. 

The input format for defining a random variable to be used as an input to the UPERT user routine 
is 

UPERT 

Notice that a DataBlock is not included, since this information will be defined in the coding for the 
UPERT subroutine. Although the statistics for the random variables may not (at present) be taken 
to be a function of time, the functional dependence of the loading on these random variables may 
be expressed as a function of time. 
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The input format for defining a perturbation variable associated with the initial velocity for a 

dynamic problem is 

VELOCITY 

nodel dvel 11 dvel21 ... dvelml 
node2 dvel 12 dvel 22 ... dvelm2 
noden dvel In dvel2n ... dvelmn 

where the values dvel are relative perturbations to the m components of the initial velocity at each 
node.. 

A slightly more cumbersome format is used to define perturbation variables associated with the 

yield stress and work-hardening curves for elastoplastic materials. However, this somewhat 

complicated format also provides a great deal of flexibility in the way elastoplastic behavior is 

defined paramaterically at different mesh nodes. The input format for defining perturbation 

variables involving the work-hardening curves is 

YIELDFUNCTION 

inodel jnodel 

dsigmal depspl 1 domegl 

dsigma2 depspl 2 domeg2 

dsigmak depsplk domegk 
inode2 jnode2 

dsigmak+1 depsplk+l domegk+1 
dsigmak+2 depsplk+2 domegk+2 

dsigma2k depspl2k domeg2k 

where dsigma, depspl, and domeg denote the relative perturbations to the yield stress, equivalent 
plastic strain, and backstress for the work-hardening curves from inode through jnode. These 
values should be arranged as a set of 3 x k tables, where k=- nhard as defined by the 
♦HARDENING option in the parameter data block. One or more node lists, each one with the 
associated work-hardening table, may be input in a single *DEFINE option, provided that all the 
tables contain exactly 3 x k entries. This concludes the discussion of the many different types of 
Level 2 perturbation variables currently available in the NESSUS finite element code. One or more 
of theses variables may be changed in a given perturbation, as discussed under the keyword 
♦PERTURB. 
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Example problem : Simply Supported Beam 


5,000 lbs 



Problem Definition: Random Variables in NESSUS input file 








PMAX 

Point Load 

loading 

5,000 lbs 

500 

Normal 

L 

Beam Length 

Geometry 

20 ft. 

2.0 

Normal 

H 

Beam Height 

Geometry 

2 ft. 

.2 

Normal 


Using the probabilistic finite element analysis program,NESSUS, determine the mean 
maximum bending stress, G x . Compare to the expected Mechanics of Solids deterministic 
solution using the flexure formula, CT=-My/I. 
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Divide your beam into 8 elements (15 nodes) with the following numbering scheme. 
Then further divide your beam to increase the accuracy of your results. 
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no * no 


♦PFEM 

C THREE POINT BEND PROBLEM 
C 

£*****************+*******+*+******************************************* 

c 

C Computational Methods for Probabilistic Engineering Analysis 
C 

Q*********************************************************************** 

c 

C NESSUS Advanced Mean Value First Order Iteration CDF analysis (AMV+) 

C 

C The material is STEEL. 

C 

C Random variables include: 

C 

C 1 PMAX point load acting at the center of the beam 

C 2 L length of the beam 

C 3 H height of the beam 

C 

C 

c 

C Z- FUNCTION DEFINITION 
C 

C 3 RANDOM VARIABLES ( 3 -COMPUTATIONAL) 

C 

♦ZFDEFINE The beginning of the z-function definition 

♦COMPUTATIONALMETHOD 1 3 This invokes the NESSUS/FEM section with 3 random 

variables . 

123 This is a list of the random variable numbers. 

C 

*END This is used to indicate the end of the ZFDEFINE section. 

C 

C 

C RANDOM VARIABLES 

C 

C 

♦RVDEFINE This is the start of the definition of the random variables. 

♦DEFINE 1 The random variable numbered as 1. 

PMAX The random variable name . 

-5000 500 NORMAL The mean, standard deviation an distribution type. 


FORCES This denotes the type of random variable. 

13 2 1.0 

DEFINE 2 This random variable is numbered 2. 

L The name of the random variable. 

20.0 2.0 NORMAL The mean, standard deviation and distribution type. 

COORDINATES The type of data block. 

1 0.00 0.0 The start of the data block. 1 is the node number 0.00 

is the perturbation amount in the x direction, 0.0 is 
the perturbation in the y direction. 

2 0.25 0.0 

3 0.50 0.0 

4 0.75 0.0 

5 1.00 0.0 

6 0.00 0.0 

7 0.25 0.0 

8 0.50 0.0 
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nnoo*o*o*on+ no 


9 

0.75 

0.0 



10 

1.00 

0.0 



11 

0.00 

0.0 



12 

0.25 

0.0 



13 

0.50 

0.0 



14 

0.75 

0.0 



15 

1.00 

0.0 



DEFINE 3 


This random variable is numbered 3. 

H 



The 

random variable name. 

2.0 

0.2 

NORMAL 

The 

mean, standard deviation and distribution type. 

COORDINATES 

The 

type of data block. 

1 

0.0 

o 

o 

The 

start of the data block. 1 is the node number 0.00 





is the perturbation amount in the x direction, 0.0 is 
the perturbation in the y direction. 

2 

0.0 

0.0 



3 

0.0 

0.0 



4 

0.0 

0.0 



5 

0.0 

0.0 



6 

0.0 

0.5 



7 

0.0 

0.5 



8 

0.0 

0.5 



9 

0.0 

0.5 



10 

0.0 

0.5 



11 

0.0 

1.0 



12 

0.0 

1.0 



13 

0.0 

1.0 



14 

0.0 

1.0 



15 

0.0 

1.0 



PERT 

1 


The 

first perturbation for random variable number one. 

1 


0.1 

The 

perturbation for the first random variable and 0.1 is the 




amount of that random variable's standard deviation that is to 




be 

perturbed. 

PERT 

2 


The 

second perturbation for random variable two. 

2 


0.1 

The 

perturbation for the second random variable and the amount 




of 

standard deviation that is to be perturbed. 

PERT 

3 


The 

third perturbation for random variable number three. 

3 


0.1 

The 

perturbation for the third random variable and the amount 




of 

standard deviation that is to be perturbed. 

END 

The 

end of 

the 

random variable parameter input. 


MEAN VALUE PROBABILISTIC ANALYSIS 

♦MVDEFINE Signals the start of the mean value analysis section. 

♦PERT 3 The number of perturbations to be probabilistically analyzed. 

123 The order of perturbations to be perturbed. 

* RAN VAR 3 The number of random variables to be used in the perturbations. 

123 The order of random variables to be perturbed. 

* DATATYPE 0 The incremental ( 0) type of data on which to perform the 
C probabilistic analysis. 

♦RESPTYPE 3 The stress type (3) of response variable to extract from the 

C perturbation database. 

♦CONDITION 0 The beginning and ending datatypes. 

♦NODE 3 The beginning and ending node numbers to probabilistically 

C analyze. Here it is the same node. 

♦COMPONENT 1 The beginnning and ending component numbers to be 

C probabilistically analyze. Here it is the x direction. 

♦PRINT LONG The keyword to signify a long printout. 
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non 


♦END The end of the MVDEFINE section 

C 

C 

C ADVANCED MEAN VALUE PROBABILISTIC ANALYSIS 
C 

♦AMVDEFINE Signals the start of the Advanced Mean Value anlysis section. 

♦ITERATION The convergence criteria for the AMV+ anlaysis. 

10 0.05 The maximum number of iterations (10) and the relative error 

C (0. 05) between consecutive runs. 

♦CONDITION 0 The beginning and ending nodes for the AMV analysis. 

♦NODE 3 The beginning and ending node numbers for the AMV analysis. 

♦COMPONENT 1 The beginnning and ending component numbers for the AMV 

C analysis. Here it is the x direction. 

♦END The end of the AMVDEFINE section. 

END PFEM INPUT 


The end of the Probabilistic Finite Element Analysis input. 


The start of the input for the deterministic finite element 
program 

This invokes the linear elastic constitutive model used for 
analysis 

This line reserves memory for 8 elements. 

This line indicates the element type is a plane stress 
element. 

This line reserves memory for 15 nodes. 

This line reserves memory for 3 boundary conditions. 

This line reserves memory for 1 applied force. 

This line prints out 10 output quantities to the ouput file. 

This line prints out 2 monitored values reported at the end of 

each iteration. 

The end of the memory allocation or parameter input. 

The iteration tolerance for the iterative solution for 3 
iterations . 

The maximum number of iterations and the relative error 
between successive iterations. 

To monitor certain values at the end of an iteration. 


STRESS 


♦COORDINATES 


TOTALDIS PLACEMENT NODE 3 COMPONENT 2 Monitor the total displacement for 

node 3 in the y direction. 

NODE 3 COMPONENT 1 Monitor the stress for node 3 in the x 

direction. 

Nodal coordinate locations 


NODE 3 COMPONENT 1 


1 

0.0 

0.0 

2 

5.0 

0.0 

3 

10.0 

0.0 

4 

15.0 

0.0 

5 

20.0 

0.0 

6 

0.0 

1.0 

7 

5.0 

1.0 

8 

10.0 

1.0 

9 

15.0 

1.0 

10 

20.0 

1.0 

11 

0.0 

2.0 

12 

5.0 

2.0 

13 

10.0 

2.0 

14 

15.0 

2.0 

15 

20.0 

2.0 


Node 1 is located at x = 0.0, y * 0.0. 
Node 2 is located at x = 5.0, y = 0.0 
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♦ELEMENTS 

151 

The element connectivity for the plane stress element 151. 

1 1 

2 7 

6 

Element number 1, nodal counter-clock-wise direction for A, 

C 



C and D. 

2 2 

3 8 

7 


3 3 

4 9 

8 


4 4 

5 10 

9 


5 6 

7 12 

11 


6 7 

8 13 

12 


7 8 

9 14 

13 


8 9 

10 15 

14 


♦BOUNDARY 


Boundary conditions 

1 2 

0.0 


Boundary at node 1 is constrained not to move in the y 

C 



direction 

1 1 

0.0 



5 2 

0.0 



♦PROPERTY 


The physical properties of the elements prescribed at the 

C 



nodes . 


1 15 1 
C 
C 
C 

* FORCES 
13 2 -5000 
* PRINT 
TOTAL NODE 
STRESS NODE 
♦END 


0 4.320E9 0.3 1.0 1.0 


The first node, last node, dummy variable, 
modulus of elasticity, Poission's ratio, 
thermal expansion coefficient and mass 
density . 

The force applied. 

Its applied at node 13 in the y direction of -5000 units. 

What quantities to print to the output file. 

Print the total displacement at all of the nodes. 

Print the stress at all of the nodes. 

The end of the model data input. 


C 

C FPI ANALYSIS CONTROL CARDS 
C 

Signifies the start of the Fast Probability Input section. 


♦FPI 

C 

♦RVNUM 3 

♦GFUNCTION 

C 

♦DATASETNM 

C 

♦METHOD 1 
C 

♦PRINTOPT i 
♦ANALTYPE ' 
♦END 
♦END 


This defines the number of random variables/perturbations. 
This defines the the Linear (1) type of g-function 
approximation . 

This defines the 4 datasets in the problem. This is necessary 
datasets are output. 

the Advanced first order method (1) solution 


as RVNUM + 1 
This defines 
technique . 
This denotes 
This defines 


the short printout (0) option. 

the full CDF analysis output (0) with 10 levels. 
This signifies the end of the FPI section. 

This signifies the end of the input file. 
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NESSUS Example problem 


The system to be analyzed is the cantilevered beam shown in Fig. 2. The beam is 
composed of steel (ASTM-A36). The system has two concentrated loads acting on it. 
Determine the maximum average normal stress and its location using NESSUS. 
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?FEM 

SFDEFUIE 

rOriPUTATIOW/VLUETIlOD 1 4 
2 3 4 
SND 


*?VDEF l MR 

DEFINE 1 

X 

00 0 . 6666G66 LOGNORMAL 
ORCES 
1 1.0 

DEFINE 2 
V 

000 0.0 LOGNORMAL 
ORCES 

2 1.0 

DEFINE 3 

.0 0.2 NORMAL 
OORD1 NATES 
0.0 0.0 
0.0 0.0 
0.0 0.0 
0.0 0.5 
. 0.0 0.5 
; o.o o.5 
• 0.0 1.0 
! 0.0 1.0 
I 0.0 1.0 

•DEFINE 4 

J 

10 0.8 NORMAL 
rOORDINATES 

1 0.0 0.0 

2 0.5 0.0 

3 1.0 0.0 

4 0.0 0.0 

5 0.5 0.0 

6 1.0 0.0 

7 0.0 0.0 

8 0.5 0.0 

9 1.0 0.0 

* PERT 1 

1 0.1 

* PERT 2 
2 0.1 

* PERT 3 

3 0.1 
•PERT 4 

4 0.1 
♦END 

C 

C 

M1VDEF1NE 
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n n 


♦ PERT 4 
12 3 4 

♦ P. ANVAR 4 
J 2 3 4 
♦DATATYPE 0 

♦ PESri YPE 3 
♦CONDITION U 
' NODE 1 

♦component l 

♦ PRINT LONG 
♦END 

C 

C 

c 

♦ AMVDEFINE 

♦ ITERATION 
10 0.05 
♦CONDITION 0 
♦NODE 1 
♦COMPONENT 1 
♦END 

C 

C 

c 

C END PFEM INPUT 
♦END 


♦FEM 
♦CONS 0 

Z . . . .MIXED METHOD 
♦ELEMENTS 4 
151 

♦NODES 9 
♦BOUNDARY 6 
♦FORCES 2 
♦PRINT 
♦MONITOR 2 

♦ END 

♦ITER 0 4 
50 0.001 
♦MONITOR 

STPESS NODE 

♦COORDINATES 

1 0.0 0.0 

2 15.0 0.0 

3 30.0 0.0 

4 0.0 1.5 

5 15.0 1.5 

6 30.0 1.5 

7 0.0 3.0 

8 15.0 3.0 

9 30.0 3.0 

* ELEMENTS 151 
112 5 4 

2 2 3 8 5 

3 4 5 8 7 

4 5 6 9 8 
♦BOUNDARY 
1 1 0.0 


COMPONENT 1 
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n n 


1 2 0.0 
4 1 0.0 
4 2 0.0 
7 1 0.0 
7 2 0.0 
» PPOPF.PTY 

1 9 J .0 J0.0E6 0.3 1.0 1.0 

* FOP'.'KS 
0 1 500 
6 2 1000 
•PR HIT 
TOTAL iiode 
STRESS MODE 

* EIJD 


*FPI 

C pfemcM 

* R VI HIM 4 
•OFUUCTIOM 1 
•DATASETMM 5 
•METHOD 1 

* PR1IITOPT 0 

* AIIALTYPE 2 
•EI0) 

•PLEVELS 11 
0.0001 
0.001 
0.01 
0.1 
* 0.2 
0.5 
0.8 
0.9 
0.99 
0.999 
0.9999 
C 

•EtlD 
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APPENDIX III 
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Executive Summary 


This report describes the work performed on the project entitled “NESSUS Enhancements and 
Technology Support.” This project was in support of the NASA grant entitled “Research and 
Education in Probabilistic Structural Analysis and Reliability” awarded to the University of 
Texas at San Antonio (UTSA). The project included NESSUS software enhancements and 
technology support for staff and students at UTSA. 

The NESSUS probabilistic structural analysis computer program combines state-of-the-art 
probabilistic algorithms with general-purpose structural analysis methods to compute the 
probabilistic response and the reliability of engineering structures. Uncertainty in loading, 
material properties, geometry, boundary conditions and initial conditions can be simulated. The 
structural analysis methods include nonlinear finite element methods, boundary element 
methods, and user-written subroutines. Several probabilistic algorithms are available such as the 
advanced mean value method and the adaptive importance sampling method. 

The previous version of NESSUS (version 6.2) was able to perform a reliability analysis of a 
structure when the failure mode is strength exceeding a stress. However, it required the user to 
develop a model for the strength portion of the limit-state function using a user subroutine. 
Another option was to define the strength (i.e., yield stress) as a single random variable. A 
material strength degradation (MSD) model has been implemented in NESSUS to provide for an 
alternative material strength capability. 

The implementation of the material strength degradation model in NESSUS was compared to 
previously published results using the PROMISS code. The PROMISS computer program 
determines the random strength of an aerospace propulsion material due to a number of random 
effects such as high and low cycle fatigue, temperature, creep, and thermal fatigue. However, 
the NESSUS results did not match the solution using PROMISS. Upon further investigation of 
NESSUS and PROMISS, it was found that the sampling scheme in PROMISS was incorrect. 
This was corrected in the updated version of PROMISS and test cases were developed and 
exercised for normal, lognormal, Weibull, and nominal variables. After this correction, 
NESSUS and PROMISS provided the same solution. 

The convergence criteria have been improved for the most probable point based methods (MPP) 
in NESSUS. Specifically for the advanced mean value method (AMV+) and the first order 
reliability method (FORM). In addition, two new keywords have been added to allow the user 
more control over the tolerances used to determine convergence. 

SwRI supported UTSA in the use and understanding of the NESSUS program. In particular, 
many questions were answered about the NESSUS input and output in support of Mark Jurena’s 
master’s thesis research. Also, questions were answered and assistance provided in developing a 
suitable problem for an undergraduate course. In addition, David Riha provided two lectures for 
the UTSA course ME 4653 “Finite Element Applications in Solid Mechanics and Design” 
describing the NESSUS computer program capabilities and several example problems (Summer 
1998 and 1999). SwRI also assisted with computer laboratories associated with these classes 
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where the students used NESSUS to compute, the probabilistic distribution of stress of a simply 
supported beam. 

Southwest Research Institute met the objectives of the grant within budget by enhancing the 
NESSUS computer program and supporting UTSA in advancing probabilistic mechanics in its 
curriculum. The NESSUS version 6.3 computer program and updated documentation were 
delivered at the completion of the project. 
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1. Introduction 


This report describes the work performed on the project entitled “NESSUS Enhancements and 
Technology Support.” This project was in support of the NASA grant entitled “Research and 
Education in Probabilistic Structural Analysis and Reliability” awarded to the University of 
Texas at San Antonio (UTSA). The scope of the project included NESSUS software 
enhancements and technology support for staff and students at UTSA. This report describes the 
accomplishments of this project 

NESSUS 


The NESSUS probabilistic structural analysis computer program combines state-of-the-art 
probabilistic algorithms with general-purpose structural analysis methods to compute the 
probabilistic response and the reliability of engineering structures. Uncertainty in loading, 
material properties, geometry, boundary conditions and initial conditions can be simulated. The 
structural analysis methods include nonlinear finite element methods, boundary element 
methods, and user-written subroutines. Several probabilistic algorithms are available such as the 
advanced mean value method and the adaptive importance sampling method. 

The application of the code includes probabilistic structural response, component and system 
reliability and risk analysis of structures considering cost of failure. Figure 1 summarizes the 
overall capabilities of NESSUS. As seen in the figure, NESSUS contains an integration of 
probabilistic methods with nonlinear finite element and boundary element methods. A general 
interface for defining random variables is included. A variety of probabilistic results can be 
obtained from the analysis of a user-defined structural model. 

Combined Stress and Strength Models 

The previous version of NESSUS (version 6.2) was able to perform a reliability analysis of a 
structure when the failure mode is strength exceeding a stress. However, it required the user to 
develop a model for the strength portion of the g-function using a user subroutine. Another 
option was to define the strength (i.e., yield stress) as a single random variable. A material 
strength degradation (MSD) model has been implemented in NESSUS to provide for an 
alternative material strength capability. All input for the MSD capability is included in the 
NESSUS/PFEM input file. All previous capabilities in NESSUS for defining the response 
function are still available. 
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Random Variables 
Loads 
-Forces 
-Pressures 
-Temperatures 

- Vtorations (PSD) 
Material Properties 

-Moduli 
-Poisson's ratio 
-Yield stress 
-Hardening parameters 

- Material orientation 
Geometry 
User-defined 

Probabilistic Methods 
Fast Probability Analysis 
-Advanced mean-value 
-First and second order 

- Fast convolution 
Sampling 

-Standard Monte Carlo 
-Latin hypercube 
-Adaptive importance 
Probabilistic Fault Tree 



Probabilistic Results 

-Full probability distribution 
-Component/single-mode reliability 
- SystenVmuttipie-modes reliability 
-Probabilistic sensItivHiea 
-ProbabiMty-based costs 

Performance Functions 
-Structural responses 

stress, strain, dlsp., freq., etc. 
-Fatigue and fracture life 
-Creep rupture Nfe 
-User-defined subroutines 
-External analysis programs 

(requires custom-made interlace) 


Analysis Types 

Static 

Transient dynamics 
Buckling 
Vibrations 
Nonlinearities 
-Plasticity 

-Large displacements 

Element Library 

Beam 

Plate 

Plane strain 
Plane stress 
Axisymmetric 
3D solid 

Enhanced solids 


Figure 1. NESSUS Capabilities 


2. Objectives and Accomplishments 

This report describes the work performed on the project entitled “NESSUS Enhancements and 
Technology Support.” This project is in support of the NASA grant entitled “Research and 
Education in Probabilistic Structural Analysis and Reliability” awarded to the University of 
Texas at San Antonio. Southwest Research Institute’s (SwRI) tasks under this project were: 


1. Provide consulting support for NESSUS-related activities, 

2. Ensure that research and educational results can be integrated into NESSUS and that 
the new capabilities are compatible with previous ones, 

3. Implement an identified suitable material model in NESSUS, 

4. Implement convergence checks into NESSUS for most probable point (MPP) based 
methods, and 

5. Correct the sampling bias in the PROMISS computer software. 

Under this contract, SwRI has successfully implemented the material strength degradation model 
and convergence checks in NESSUS and the associated documentation has been updated for the 
new capabilities. The PROMISS computer code was modified to correct a sampling bias. 
Consulting support for NESSUS-related activities was provided to UTSA researchers and SwRI 
assisted UTSA in instructing two undergraduate courses dealing with NESSUS. These tasks 
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were completed within the allocated funding. The following sections describe the completed 
tasks. 

3. Material Model Implementation in NESSUS 

3.1 Material Strength Degradation Model 

The multi-factor equation for material strength degradation (MSD) has been implemented in 
NESSUS using the NZFUNC subroutine that was initially developed for predefined resistance 
models. All input for the MSD model is in the NESSUS/PFEM input deck using a keyword 
interface consistent with previous versions of NESSUS. The required input is fully described in 
an updated version of the NESSUS/PFEM User’s manual. In addition, all previous capabilities 
for defining the response function in NESSUS are still available. 

This MSD model is based on the multi-factor equation defined in Reference 2. This is an 
empirical equation and the coefficients are defined using test data and regression. Reference 2 
provides values for Inconel 718. The model is currently set up for five fixed effects and 12 user 
defined effects. These effects are 

• Temperature 

• High cycle mechanical fatigue 

• Low cycle mechanical fatigue 

• Creep 

• Thermal fatigue 

• Up to 12 user defined effects 

The form of the equation is 

'j^rz£ t r h zi t \Kzkt r 4 ,-a t 

L t v -to] Ix-tfoj Ia-V 

Temperature HCF LCF Creep Thermal Fatigue User 

Where T is temperature, N is cycles, t is time, and A is used for user defined generic quantities. 
The lower case exponents are the empirical material constants for each effect. The U subscript 
refers to an ultimate value, the O subscript is the reference value, and the unscripted term is the 
current value for the effect. 

To increase model sensitivity for any effect as discussed in Reference 2, a log transformation can 
be introduced. For Inconel 718, all effects except temperature use the log transformation and 
have the form: 


5 =5 Li Zi \IjlzJL 

°T U -T 0 Nfj - N 0 
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s = s 0 


LOG (Aw ) - LOG (A f ) T 

LOG ( A, y ) - LOG (A,o ) J 

\ 


For other materials, the nature of the data will dictate the use of the log transformation. The 
effects to be used and the model type (log transformation) are defined in the Z-function 
definition section of the PFEM input file. Any combination of effects can be selected and any of 
the terms can be considered random if desired. 


The random variables used in MSD model can also be used in the computational model. For 
example, if temperature is a random variable, it can be a temperature loading on the finite 
element model (contributing uncertainty to the stress) and a term in the material strength model 
(contributing uncertainty to the strength). 

The form of the Z-function used by NESSUS is 


z = 



- a 


where 

So is the reference value of material strength 

Am is the ultimate value of the particular effect 

Aio is the reference value of the particular effect 

Ai is the current value of the particular effect 

a, is the empirical material constant for the particular effect 

a is the structural response, i.e., stress 


3.2 Material Strength Degradation Model Input 

The material strength degradation model input is described in the NESSUS/PFEM user’s 
manual. The section on predefined resistance models describes the MSD model and all keywords 
used are defined in the keyword section of the manual. An example problem is presented in the 
next section to assist with using this new capability. 


33 Example Problem 

The example problem is to compute S/So for two effects, high cycle mechanical fatigue and high 
temperature. The equation for these two effects is 


s 

i 

i 

£ 
i 

9 

' LOG (N„) - LOG 

■nw 

(N) 

s 0 ~ 

1 

£ 

1 

£ 
i 


LOG (N v ) - LOG 

(tfo). 
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The formulation defined in NESSUS is 


Z = S 0 




LOG (N 0 ) ~ LOG (AQ 
LOG (N v ) - LOG ( N 0 ) 




- a 


By defining 5 to be 1.0 and not including a computational model to compute stress (a = 0.0), 
then NESSUS is solving the problem for S/S 0 . The random variables are defined in Table 1 


Table 1. Random Variable Definitions 


EfTect 


Distribution 

Mean 

Standard 

Deviation 

High Cycle 

Nu 

Normal 

1. Ox 10 10 cycles 

1. Ox 10 9 cycles 

Mechanical 

N 

Normal 

2.5x10 s cycles 

2.5x1 0 4 cycles 

Fatigue 

No 

Normal 

0.25 cycles 

0.025 cycles 

(at 75 °F) 

W 

Normal 

0.3785 

0.0114 

High 

Tu 

Normal 

2369.0° 

236.90° 

Temperature 

T 

Normal 

1000.0° 

100.0° 

(at 1000 °F) 

To 

Normal 

75.0° 

7.5° 


0 

Normal 

0.2422 

0.0088 


The relevant NESSUS input for this problem is shown in Figure 2. Note that to consider a term 
in an effect random, the alphanumeric name must match the name in the RVDEFINE input. An 
example of this is the highlighted NU variable shown in both the ZFDEFINE section and the 
RVDEFINE section in Figure 2. 
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•ZFDEFINE 
•EXPLICITMODEL 8 
1.2. 3. 4, 5, 6,7, 8 
•ZFUNCT 3 0 

•MSDM 1.0 (So) 

TEMP LINEAR TU T TO QQ 
HCF LOG MO N NO S 
END 

*END 

•RVDEFINE 
•DEFINE 1 

MO 

1.0E10 1.0E9 NORMAL 
•DEFINE 2 
N 

2.5E5 2.5E4 NORMAL 

•DEFINE 3 
NO 

0.25 0.025 NORMAL 

•DEFINE 4 
S 

0.3785 0.0114 NORMAL 

•DEFINE 5 
TU 

2369.0 236.9 NORMAL 
•DEFINE 6 
T 

1000.0 100.0 NORMAL 
•DEFINE 7 
TO 

75.0 7.50 NORMAL 

•DEFINE 8 
QQ 

0.2422 0.0088 NORMAL 

Figure 2. NESSUS Input for Example Problem 


The cumulative distribution function is shown in Figure 3. 
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Lifetime Strength 


— NESSUS AMV+ ♦ PROMISS 



0.45 0.5 0.55 0.6 0.65 0.7 0.75 


Lifetime Strength 

— NESSUS AMV+ ♦ PROMISS 


Figure 3. CDF for Example Problem 


4. Convergence Checks Implementation in NESSUS 

The convergence criteria have been improved for the most probable point based methods (MPP) 
in NESSUS, specifically for the advanced mean value method (AMV+) and the first order 
reliability method (FORM). In addition, two new keywords have been added to allow the user 
more control over the tolerances used to determine convergence. Default values for all 
tolerances are still available but there are times when these may need to be adjusted. The first 
order reliability method (FORM) in NESSUS is the cornerstone of several other methods. The 
following methods use the most probable point (MPP) as a starting point for the probability 
calculation: 

second order reliability method (FORM) 

fast probability integration (FPI) 

linear based adaptive importance sampling (AIS1) 

curvature based adaptive importance sampling (AIS2) 

Harbitz importance sampling (ISAMF) 

The FORM or FPI method is also used in the AMV+ procedure to compute the probability based 
on the approximate function. 

The FORM uses the Rackitz-Feissler optimization algorithm to locate the MPP. While this is an 
efficient method and works well for most well behaved problems, it is not guaranteed to 
converge. There have been several instances where the FORM solution has not calculated the 
correct results. In many cases this is because of an unobtainable performance value. The 
previous version of NESSUS did not check the Z value to insure that the located MPP was on the 
limit-state function. The p-level algorithm using FORM as implemented in NESSUS uses a 
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quadratic fit of the CDF to determine the. Z value for the input probability value. New 
convergence checks and additional control over the convergence checks have been implemented. 

Similar convergence criteria have been implemented in NESSUS for the AMV+ and FORM 
algorithms. The convergence criteria is different when using the p-level (user inputs a probability 
level and NESSUS computes the performance value) or z-level (user inputs a performance value 
and NESSUS computes the probability value) algorithms. Further details about the convergence 
criteria and new keyword options are contained in the updated NESSUS documentation. The 
following convergence criteria are implemented in NESSUS for the AMV+ and FORM 
algorithms. 

FORM Convergence Criteria 

For the z-level algorithm convergence is defined as 

a) the relative change in P is less than a default (0.0001) or user specified tolerance, 
beta_tol. 


A 

A-i 


^ beta_tol 


(note that beta_tol is used for P<4.0 and 10* beta_tol for P>4.0) 


AND 


b) the difference in the computed value of Z and the input value of Z* is less than a 
maximum allowable tolerance, ztol, times the approximate standard deviation, o ^ 
(the default is 0.001) 


AND 


Z-Z* <ztol-a w 


c) the measure of the angle between the MPFs from the last two iterations, 6, is less 
than a maximum allowable (default of 30° or user defined). 6 is defined as cos 0 = oti 
• 0Cj.i where oti are the direction cosines to the MPP from iteration i. 

For the p-level algorithm convergence is defined as 

a) the relative change in z is less than a default (.001) or user specified tolerance, ztol, 
times the approximate standard deviation , a , 


AND 


AZ <ztol*<r ¥p 
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b) the measure of the angle between the MPP's from the last two iterations, 0, is less 
than a maximum allowable (default of 30° or user defined). 0 is defined as cos 0 = Gj 
• otj.i where a* are the direction cosines to the MPP from iteration i. 


AMV+ Convergence Criteria 

For the z-level algorithm convergence is defined as 

a) the relative change in [3 is less than a user specified tolerance, beta_tol (default is 

0 . 01 ) 


beta_tol 


AND 


b) the relative change in the computed value of Z compared to the input value of Z* is 
less than a maximum allowable, ztol (default is 0.01), 


AND 



< z tol 


c) the measure of the angle between the MPP's from the last two iterations, 0, is less 
than a maximum allowable (default of 30° or user defined). 0 is defined as cos 0 = otj 
• 0,-1 where cti are the direction cosines to the MPP from iteration i. 


For the p-level algorithm convergence is defined as 

a) the relative change in z is less than a default (0.01) or user specified tolerance, z_tol. 


AND 



<, z_tol 


b) the measure of the angle between the MPP's from the last two iterations, 0, is less 
than a maximum allowable (default of 30° or user defined). 0 is defined as cos 0 = otj 
• Oj.i where ou are the direction cosines to the MPP from iteration i. 
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5. PROMISS Code Update 


The PROMISS computer program determines the random strength of an aerospace propulsion 
material due to a number of random effects such as high and low cycle fatigue, temperature, 
creep, and thermal fatigue. The implementation of the material strength degradation model in 
NESSUS was verified using data for Inconel 718 (Reference 2) and compared to the cumulative 
distribution function. The NESSUS computed CDF did not match the results using PROMISS. 
Upon further investigation of NESSUS and PROMISS, it was found that the sampling scheme in 
PROMISS was incorrect. The same starting seed was used for sampling each random variable. 
This was corrected in the updated version of PROMISS and test cases were developed and 
exercised for normal, lognormal, Weibull, and nominal variables. In addition for testing 
purposes, input and output files were given a common file name with appropriate extensions. 

6. Support of NESSUS Activities at UTSA 

Southwest Research Institute (SwRI) supported UTSA in the use and understanding of the 
NESSUS program. In particular, many questions were answered about the NESSUS input and 
output in support of Mark Jurena’s master’s thesis research. Also, questions were answered and 
assistance provided in developing a suitable problem for an undergraduate course. 

7. Support of Educational Activities at UTSA 

Southwest Research Institute (David Riha) provided two lectures for the UTSA course ME 4653 
“Finite Element Applications in Solid Mechanics and Design” describing the NESSUS computer 
program capabilities and several example problems (Summer 1998 and 1999). SwRI also 
assisted with computer laboratories associated with these classes where the students used 
NESSUS to compute the probabilistic distribution of stress of a simply supported beam. 

8. Conclusions 

Previous versions of NESSUS required the user to develop a user subroutine to model the 
strength portion for a structural reliability analysis or use a single random variable for the 
strength (i.e., yield stress). The material strength degradation model (MSD) has been 
successfully implemented and tested in the NESSUS computer program. The MSD model 
provides a new capability in NESSUS in that there is now a material strength model option for 
performing reliability analysis with NESSUS that is defined by keywords in the NESSUS input 
file. All previous capabilities for defining the response function are still available in NESSUS. 
After correcting a sampling bias in the PROMISS code, NESSUS and PROMISS provide the 
same results. 

The convergence criteria have been improved for the most probable point based methods (MPP) 
in NESSUS. Specifically for the advanced mean value method (AMV+) and the first order 
reliability method (FORM). In addition, two new keywords have been added to allow the user 
more control over the tolerances used to determine convergence. These improved convergence 
checks will help prevent NESSUS from providing wrong results. This information could be used 
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to switch to another method to locate the MPP. Future enhancements to NESSUS may include 
automatically switching to another MPP search method when the solution does not convergence. 

Southwest Research Institute also aided the researchers at UTSA with the use of NESSUS in 
explaining the input and output of the NESSUS program. In addition, SwRI assisted with two 
undergraduate courses at UTSA by providing lectures about NESSUS and assisting with 
computer laboratories in which the NESSUS program was used. 
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1 Introduction 


These release notes describe the modifications and enhancements of version 6.3 of the NESSUS 
probabilistic structural analysis code. The relevant sections of the NESSUS documentation have 
been updated for the new features. The new features in version 6.3 include: 

• Material Strength Degradation Model 

• MPP (AMV+ and FORM) Convergence Enhancements 

• Other Modifications and Enhancements 

The following sections describe these enhancements and modifications in more detail. The 
appendix at the end of this document lists the modified NESSUS subroutines. 


2 Material Strength Degradation Model 

The previous version of NESSUS (version 6.2) was able to perform a reliability analysis of a 
structure when the failure mode is strength exceeding a stress. However, it required the user to 
develop a model for the strength portion of the limit-state function using a user subroutine. 
Another option was to define the strength (i.e., yield stress) as a single random variable. A 
material strength degradation (MSD) model has been implemented in NESSUS to provide for an 
alternative material strength capability. 

The multi-factor equation for material strength degradation (MSD) has been implemented in 
NESSUS using the NZFUNC subroutine that was initially developed for predefined resistance 
models. All input for the MSD model is in the NESSUS/PFEM input deck using a keyword 
interface consistent with previous versions of NESSUS. The required input is fully described in 
an updated version of the NESSUS/PFEM User’s manual. In addition, all previous capabilities 
for defining the response function in NESSUS are still available. 

This MSD model is based on the multi-factor equation defined in Reference 1. This is an 
empirical equation and the coefficients are defined using test data and regression. Reference 1 
provides values for Inconel 718. The model is currently set up for five fixed effects and 12 user 
defined effects. These effects are 

• Temperature 

• High cycle mechanical fatigue 

• Low cycle mechanical fatigue 

• Creep 

• Thermal fatigue 

• Up to 12 user defined effects 
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The form of the equation is 
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Where T is temperature, N is cycles, t is time, and A is used for user defined generic quantities. 
The lower case exponents are the empirical material constants for each effect. The if subscript 
refers to an ultimate value, the O subscript is the reference value, and the unscripted term is the 
current value for the effect. 

To increase model sensitivity for any effect as discussed in Reference 1 , a log transformation can 
be introduced. For Inconel 718, all effects except temperature use the log transformation and 
have the form: 

LOG (A lV ) - LOG {A, ) T 
"LlOG(4„)-LOG(^„) 

For other materials, the nature of the data will dictate the use of the log transformation. The 
effects to be used and the model type (log transformation) are defined in the Z-function 
definition section of the PFEM input file. Any combination of effects can be selected and any of 
the terms can be considered random if desired. 

The random variables used in the MSD model can also be used in the computational model. For 
example, if temperature is a random variable, it can be a temperature loading on the finite 
element model (contributing uncertainty to the stress) and a term in the material strength model 
(contributing uncertainty to the strength). 

The form of the Z-function used by NESSUS is 


z = 



- a 


where 

So is the reference value of material strength 

Aiu is the ultimate value of the particular effect 

Aio is the reference value of the particular effect 

Ai is the current value of the particular effect 

a ( is the empirical material constant for the particular effect 

a is the structural response, i.e., stress 
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2. 1 Material Strength Degradation Model Input 


The material strength degradation model input is described in the NESSUS/PFEM user’s 
manual. The section on predefined resistance models describes the MSD model and all keywords 
used are defined in the keyword section of the manual. An example problem is presented in the 
next section to assist with using this new capability. 


2.2 Example Problem 

The example problem is to compute S/So for two effects, high cycle mechanical fatigue and high 
temperature. The equation for these two effects is 


5 

1 

1 

1 

<7 

' LOG (N u ) - 

LOG (N) ' 

So" 

1 

<? 

1 

1 


LOG (Ny) - 

LOG (N a ) 


The formulation defined in NESSUS is 


Z = S r 


& 

1 

1 

1 

' LOG ( Ny ) - LOG (N) " 

- MUBm t 

1 

E-? 

LOG (Ny) - LOG (N a ) 


- a 


By defining S to be 1 .0 and not including a computational model to compute stress (a = 0.0), 
then NESSUS is solving the problem for S/So • The random variables are defined in Table 1. 

Table 1. Random Variable Definitions 


Effect 

Symbol 

Distribution 

Mean 

Standard 

Deviation 

High Cycle 

Nu 

Normal 

1. Ox 10 lu cycles 

1. Ox 10 v cycles 

Mechanical 

N 

Normal 

2.5x10 s cycles 

2.5xl0 4 cycles 

Fatigue 

No 

Normal 

0.25 cycles 

0.025 cycles 

(at 75 °F) 

W 

Normal 

0.3785 

0.0114 

High 

Tu 

Normal 

2369.0° 

236.90° 

Temperature 

T 

Normal 

1000.0° 

100.0° 

(at 1000 °F) 

To 

Normal 

75.0° 

7.5° 


Q 

Normal 

0.2422 
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The relevant NESSUS input for this problem is shown in Figure 1. Note that to consider a term 
in an effect random, the alphanumeric name must match the name in the RVDEFINE input. An 
example of this is the highlighted NU variable shown in both the ZFDEFINE section and the 
RVDEFINE section in Figure 1. 


‘ZFDEFINE 
‘EXPLICITMODEL 8 
1,2, 3, 4, 5, 6, 7, 8 
‘ZFUNCT 3 0 

‘MS DM 1.0 (S 0 ) 

TEMP LINEAR TU T TO QQ 
HCF LOG NU N NO S 
END 

‘END 

‘RVDEFINE 
‘DEFINE 1 
NU 

1.0E10 1.0E9 NORMAL 
‘DEFINE 2 
N 

2.5E5 2.5E4 NORMAL 

‘DEFINE 3 
NO 

0.25 0.025 NORMAL 

‘DEFINE 4 
S 

0.3785 0.0114 NORMAL 

‘DEFINE 5 
TU 

2369.0 236.9 NORMAL 
‘DEFINE 6 

T 

1000.0 100.0 NORMAL 
‘DEFINE 7 

TO 

75.0 7.50 NORMAL 

‘DEFINE 8 

QQ 

0.2422 0.0088 NORMAL 


Figure 1. NESSUS Input for Example Problem 


The cumulative distribution function is shown in Figure 2. 
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Figure 2. CDF for 
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Example Problem 


3 MPP (AMV+ and FORM) Convergence Enhancements 

The convergence criteria have been improved for the most probable point based methods (MPP) 
in NESSUS, specifically for the advanced mean value method (AMV+) and the first order 
reliability method (FORM). In addition, two new keywords have been added to allow the user 
more control over the tolerances used to determine convergence. Default values for all 
tolerances are still available but there are times when these may need to be adjusted. The first 
order reliability method (FORM) in NESSUS is the cornerstone of several other methods. The 
following methods use the most probable point (MPP) as a starting point for the probability 
calculation: 

second order reliability method (FORM) 

fast probability integration (FPI) 

linear based adaptive importance sampling (AIS1) 

curvature based adaptive importance sampling (AIS2) 

Harbitz importance sampling (ISAMF) 

The FORM or FPI method is also used in the AMV+ procedure to compute the probability based 
on the approximate function. 

The FORM uses the Rackitz-Feissler optimization algorithm to locate the MPP. While this is an 
efficient method and works well for most well behaved problems, it is not guaranteed to 
converge. There have been several instances where the FORM solution has not calculated the 
correct results. In many cases this is because of an unobtainable performance value. The 
previous version of NESSUS did not check the Z value to insure that the located MPP was on the 
limit-state function. The p-level algorithm using FORM as implemented in NESSUS uses a 


N AS A/CR— 2001-211112 


121 





quadratic fit of the CDF to determine the Z value for the input probability value. New 
convergence checks and additional control over the convergence checks have been implemented. 

Similar convergence criteria have been implemented in NESSUS for the AMV+ and FORM 
algorithms. The convergence criteria is different when using the p-level (user inputs a probability 
level and NESSUS computes the performance value) or z-level (user inputs a performance value 
and NESSUS computes the probability value) algorithms. Further details about the convergence 
criteria and new keyword options are contained in the updated NESSUS documentation. The 
following convergence criteria are implemented in NESSUS for the AMV+ and FORM 
algorithms. 


3.1 FORM Convergence Criteria 

For the z-level algorithm convergence is defined as 

a) the relative change in P is less than a default (0.0001) or user specified tolerance, 
beta_tol. 


I A 

A-, 


beta_tol 


(note that betatol is used for P<4.0 and 10* betatol for p>4.0) 

AND 


b) the difference in the computed value of Z and the input value of Z* is less than a 
maximum allowable tolerance, ztol, times the approximate standard deviation a 

5 °PP 

(the default is 0.001) 


AND 


Z-Z* <ztol-cr 

w 


c) the measure of the angle between the MPP's from the last two iterations, 0, is less 
than a maximum allowable (default of 30° or user defined). 0 is defined as cos 0 = Gq 
• a;, i where oq are the direction cosines to the MPP from iteration /. 

For the /7-level algorithm convergence is defined as 

a) the relative change in z is less than a default (.001) or user specified tolerance, ztol, 
times the approximate standard deviation , , 


AND 


AZ 

VP 


NASA/CR— 2001-211112 


122 



b) the measure of tire angle between the MPP's from the last two iterations, 0, is less 
than a maximum allowable (default of 30° or User defined). 0 is defined as cos 0 = aj 
• aw where ai are the direction cosines to the MPP from iteration /. 


3.2 AMV+ Convergence Criteria 

For the z-level algorithm convergence is defined as 

a) the relative change in P is less than a user specified tolerance, beta_tol (default is 

0 . 01 ) 


I fl-ft., L 


p. 


beta tol 


AND 


b) the relative change in the computed value of Z compared to the input value of Z* is 
less than a maximum allowable, ztol (default is 0.01), 


Z-Z‘ 

^ z tol 

Z 

AND 


c) the measure of the angle between the MPP's from the last two iterations, 0, is less 
than a maximum allowable (default of 30° or user defined). 0 is defined as cos 0 = aj 
• aj.| where ai are the direction cosines to the MPP from iteration i. 

For the p-level algorithm convergence is defined as 

a) the relative change in z is less than a default (0.01) or user specified tolerance, z_tol, 


AND 



< z tol 


b) the measure of the angle between the MPP's from the last two iterations, 0, is less 
than a maximum allowable (default of 30° or user defined). 0 is defined as cos 0 = a* 
• aj.| where aj are the direction cosines to the MPP from iteration i. 
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4 Other Modifications and Enhancements 

The following lists identifies several other corrections and enhancements implemented in 
NESSUS version 6.3: 

• corrected error check on Z for AMV+ analysis 

• corrected error for implicit sampling using Harbitz method 

• changed default tolerance from 0.2 to 0.01 on Z check for AMV+ 

• tightened tolerance on chi square distribution iteration 

• skip integration for special cases of the maximum entropy distribution (e.g. uniform) 

• check for valid FPI methods for AMV analyses, correlated variables and implicit 
sampling 

• corrected open statement error in SIMFEM 


5 References 


1. Bast, C. and Boyce, L., “Probabilistic Material Strength Degradation Model for Inconel 718 
Components Subjected to High Temperature, High-Cycle and Low-Cycle Mechanical 
Fatigue, Creep and Thermal Fatigue Effects,” NASA CR 198426, NASA Lewis Research 
Center, Cleveland, Ohio, November 1995. 
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6 Appendix 

The following tables describe the modified subroutines for NESSUS 6.3. 

Table 2. Modified NESSUS Subroutines for Port from Cray to SGI 


Description 

routines 

modified 

replaced CPU time system call for SGI 

quit.f 

date and time system call on SGI 

dater.f 

double precision flag change for 32 bit machine (IDP) 

nessus.f 

modified arguments as required on SGI for call to GETARG 

prompt.f 

replaced CPU time system call for SGI 

timer.f 

initialize ICREAD=5 in nessus main routine 

nessus.f 

set NDBREC to 80 (file record size) 

in tint. f 


Table 3. Modified NESSUS Subroutines for the Material 
Strength Degradation Model (MSDM) 


Description 

routines 

modified 

added predefined equation 3 for material strength degradation model 

nzfunct.f 

added call to rdmfie to read MSDM input 

redpfm.f 

new routine to read MSD model input 

rdmfie. f 

added call to msdump 

pfdump.f 

new routine to echo MSD model input 

msdump.f 

added test model 

respon.f 


Table 4. Modified NESSUS Subroutines for AMY Convergence Enhancements 


Description 

Routines 

modified 

added input echo for *CONVERGENCE_CONTROL card 

pfdump.f 

added user control over convergence tolerances 

pfemal.f 

added user control over convergence tolerances 

pfnopt.f 

added user control over convergence tolerances 

pfyspt.f 

added input and checks *CONVERGENCE_CONTROL card 

redpfm.f 

added user control over convergence tolerances 

rfemal.f 

added user control over convergence tolerances 

uconv.f 

added write statement for new convergence tolerances 

wrtmov.f 
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Table 5. Modified NESSUS Subroutines for FORM Convergence Enhancements 


Description 

routines 

modified 

added z and 9 convergence checks and user control over tolerances 

fit.f 

added input for the *FCONVERGENCE card 

redmod.f 

added z and 9 convergence checks and user control over tolerances 

zlevel.f 

added write of *FCONVERGENCE card to FPI input deck 

fpiter.f 


Table 6. Modified NESSUS Subroutines for Other Modifications 


Description 

routines 

modified 

Changed version number and release date 

nessus.f 

Corrected format for year 2000 

header, f 
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Table 7. Modified NESSUS Subroutines for Code Corrections 


Description 

Bus error when using Harbitz method 
through PFEM. Only affects jobs 
using ISAMF and IMPLICIT during a 
PFEM run. 


Correction 

RWORK array not dimensioned 
in FASTMC. Add RWORK to all 
calls to FASTMC and other 
calling subroutines. 


Error check on Z incorrect when Z is 

negative for AMV+ analysis 

Improved tolerance 


Take absolute value of ZFLTN 

when making check. 

change ZERR_TOL form 0.2 to 


routines 

modified 

cdfglb.f 

cdfloc.f 

cintvl.f 

fastmc.f 

xfpi.f 

zlevel.f 

uconv.f 

uconv.f 


Allow small negative eigen values 
when transforming from correlated 
variables to independent variables. 
Can appear for very highly correlated 
random variables where some 
independent variables are expected to 
have variance at or very close to zero. 
Caused by the numerical algorithm. 
Tighten tolerance on chi square 
distribution 


Skip integration of maximum entropy 
distribution if the pdf is constant 

(uniform distribution). 

Checks for valid FPI methods for 
AMV analyses, correlated variables, 

and implicit sampling. 

Open too many files on some systems 

when running SIMFEM 

Passes single precision number as 
double precision to XINV (for 

confidence interval) 

Incorrectly made AD a small number 


0.01 

Set the eigen value to 10e-10 
times the maximum eigen value 
when negative and if the absolute 
value is within 1/1000 of the 
maximum eigen value value. 


change convergence criteria from 
IE-4 to IE-8 and increase 
iterations from 10 to 50. 

Add check for constant 
distribution 

Added checks after reading input. 


add a close statement for the 
scratch file 
set variable P95 


removed statement 


rottrn.f 


chi.f 


cdf6.f 


redpfm.f 


simres.f 

fastmc.f 


decomp, f 
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Abstract 

This paper deals with stress analysis of a cantilever beam. The analysis was performed 
three different ways. The first analysis of the stresses at different regions within the 
body were done using simple beam theory. This is an algebraic way of arriving at the 
stresses without having to do extensive calculations. Stress components were calculated 
at the location of maximum stress. The maximum stress within the beam was found to 
occur at the fixed end, it happened to be the normal stress in the longitudinal direction of 
the beam. Its value was calculated to be 20.2 k.s.i. The second way the stress was 
calculated was by modeling and analyzing the same beam using NESSUS, a probabilistic 
finite element code. NESSUS accounts for model uncertainty by allowing the user to 
input one or more random design variables. Random variables have a range of possible 
values and associated probabilities of occurring. Three parameters were modeled as being 
random: the beam’s length, its height, the horizontal load, and the vertical load. Thus, the 
output will also exhibit uncertainty. Three runs were performed using NESSUS. For the 
first run, the beam was discretized into 4 plane stress elements. The displacement 
method of analysis was used; hence, NESSUS assumed a displacement field form. The 
maximum stress occurred at the fixed end and was found to be the normal stress in the x- 
direction. It’s average value was 15 k.s.i, which is a 25.7% difference with respect to 
simple beam theory. For the second run, the beam was modeled with 40 plane stress 
elements. The displacement method of analysis was also used for this run. The 
maximum stress was noted to occur in the same location and was the normal stress in the 
x-direction. It’s average value was 19.5 k.s.i, which is a 3.5% difference with respect to 
simple beam theory. The third and final way of analyzing the beam in NESSUS was by 
modeling the beam with 4 plane stress elements, except the mixed method of analysis was 
used by NESSUS. This means that a stress and a displacement field form were assumed 
by the code. The maximum stress occurred in the same location and was in the same 
direction. It’s average value was 20.3 k.s.i, which is a 0.5% difference with respect to 
simple beam theory. The final way of calculating the stresses within the beam was done 
by writing a finite element code, which used the finite element method by way of 
weighted residuals on the plane stress elasticity equations. The displacement method was 
used, and the beam was discretized into 40 plane stress elements. The Galerkin-Bubov 
approach was implemented; hence, the weight functions were chosen from the same 
shape functions used to bilinearly approximate the displacement field vector. Also, 
isoparametric mapping was used; thus, because bilinear shape functions were used to 
approximate the displacement field, the global coordinates were bilinearly mapped to the 
local coordinates. The maximum stress occurred at the same location and is also the 
normal stress in the x-direction. It’s value was calculated to be 18.8 k.s.i, which is a 
3.59% difference with respect to the 40 element, displacement method run performed by 
NESSUS. A cumulative distribution function was obtained for the 3rd mixed method run. 
From this probabilistic output, one was able to deduce that there was a 99.99% 
probability that the maximum normal stress in the x-direction will be less than or equal to 
30.4 k.s.i. Therefore, there is a probability smaller than 1/10000 that the beam will fail by 
yielding. 
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Introduction 

A Cantilevered beam is useful model of reality that engineers use to analyze real 
structures or parts of real structures. Centrifugal pump impellers and Space Shuttle Main 
Engine turbopump turbine blades are examples of structures that can be initially modeled 
as cantilevered beam. Thus, a cantilevered beam is a good, simple model of reality that is 
useful to acquire a “feel” for the stress magnitudes within a structure. 

Many times one would analyze a whole or part of a system to ensure that failure 
does not occur in a way that hinders the performance of the system. If a system does not 
perform well when a piece of it breaks, then that is the way it fails. If a system does not 
perform well when a component permanently deforms or deforms too much, then that is 
its respective mode of failure. For this paper, we will ensure that a cantilevered beam 
does not fail by permanently deforming, or yielding. Thus, we analyze to ensure that a 
system does not fail in a certain way. The system to be analyzed here is a cantilevered 
beam because it is a good model of things like a turbine blade, which is expensive; and, at 
their speed of operation, some types of failure can be disastrous. 


The System 

The system that was analyzed was the cantilevered beam shown in Figure 1 . A 
characteristic of a 
cantilevered beam is that the 
horizontal and vertical 
displacements are zero for 
all points on the fixed end, 
this is known as a boundary 
condition. The beam was 
composed of steel (ASTM- 
A36). This material has 
several associated average 
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material properties. For calculations, an average modulus of elasticity and an average 
Poisson’s ratio of 30,000 k.s.i. and 0.3, respectively, were used. The materials average 
yield strengths in tension and in shear were found to be 36 and 21 k.s.i., respectively. 

The system also has certain average geometric properties: its depth is 1.0 inch, the 
height is 3.0 inches, the length is 30 inches, and its area moment of inertia about the axis 
normal to the page was calculated from the average dimensions to be 2.25 in' . 

The system has two concentrated loads acting upon it, as shown in Figure 1. 

These loads cause the system to deform and exhibit internal stresses. A horizontal load 
and a vertical load, whose average values are 500 lbs. and 1000 lbs., respectively, are 
applied to the beam’s free end, at a point halfway up the height of the beam and 
uniformly distributed across the thickness. Also, both loads are assumed to be applied 
quasi-statically so that the dynamic effects of the load application are to be ignored; and, 
the beam is considered to be in static equilibrium throughout its analysis. 


NESSUS 

The analysis for this system was performed using NASA Glenn Research 
Center’s Numerical Evaluation of Stochastic Structures Under Stress (NESSUS) software, 
a probabilistic finite element code. It has the capability of analyzing static and dynamic 
problems, linear or even nonlinear problems, and it can compute the sensitivity of the 
wanted output calculations with respect to small variations in the user defined random 
input variables. 

Three different plane stress analyses were accomplished using NESSUS. 
Therefore, the following assumptions of plane stress must be true: 

1. Only 0^,0^, and o „ are nonzero, the other 6 stress components are zero. 

2. Any body force in the z direction is assumed to be zero. 

3. As a result of these assumptions and Hooke’s Law, which relates stress to strain, we 
find that strains in the z direction ( y xz and y yz ) are also zero. 
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4. All variables, stresses, forces, strains, etc. are functions of x and y only, and any 
loading is assumed to be distributed uniformly across the thickness. 

5. Also, from it is known that plane stress prevails in a beam that is thin with respect to 
its other dimensions [3]. 

Therefore, since our beam shown in Figure 1 meets all of the assumptions previously 
stated (as well as those for simple beam theory, which will be discussed later), the plane 
stress elements in NESSUS can be used to model our beam; and, one should expect results 
close to those obtained using simple beam theory. 

NESSUS uses a finite element method to obtain a solution. It is an approximate 
method because one cannot generally solve the differential equations that will give us all 
of the stress or deformation components at every point within the body. Thus, we give 
up on the hopes of an exact solution and solve the equations only at certain points within 
the body. These points are known as nodes. Based on the values of our solution at these 
nodes, we can estimate the solution values between the nodes. 

The NESSUS code lets one enter certain parameters of our physical model as 
random variables. Random variables are variables that can be any of a set of values, be 
this set continuous (infinite in amount), or discrete; and, each value has an associated 
chance or probability of occurring. These values and their respective chances or 
probabilities of occurring form a probability distribution function (p.d.f.) for that random 
variable. Table 1 shows the probabilistic input parameters to be used for all runs in 
NESSUS. 

Table 1. Probabilistic input parameters of the cantilever beam. 


Name 

Symbol 

Distribution 
Type (mean, 
std. dev.) 

Units 

Length 

L 

Normal 
(30.0, 0.8) 

in. 

Height 

h 

Normal 
(3.0, 0.2) 

in. 

Horizontal 

P 

Lognormal 

lb. 

Load 

r X 

(500, 0.666) 


Vertical 

P 

Lognormal 

lb. 

Load 

r y 

(1000, 0.8) 
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Note the distribution types used for geometric and load variables and that these 
p.d.f.s need only the mean value and the standard deviation to completely describe the 
respective random variable. Using these random variables and the deterministic or single 
valued parameters shown in Figure 1, several different analyses were performed using 
NESSUS.' Some typical outputs that NESSUS can compute are the displacements, 
velocities, accelerations, stresses and strains for each node. For this report, the stress is 
the output variable and it will be in the form of a random variable. The stress output at 
each user requested node will have a range of values and each value will have an associated 
probability of occurring (an output p.d.f. is the result). 

1st run. 

The first run in NESSUS had a node scheme like that shown in Figure 2. The system 
was discretized into 4 elements and 9 nodes. The displacement method was used by the 
software to solve the system of equations, which means that the displacement field took 
on an assumed form. This 
analysis required a stress 
output for each node to ensure 
that the beam does not fail by 
permanently deforming. For 
this first run, the maximum 
average normal stress on the x 
face was found to be 
approximately 15.0 k.s.i., it 
occurred at node 1 . Also, all 

shear stresses within the beam were less than 1 k.s.i. 
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2nd run. 

The second run in NESSUS had the same nodal scheme as the first run. The system 

v. * w „ 

model can be seen in Figure 3. This time, the'mixed method of analysis was used by the 
software in solving the system of equations. This means that the software assumes a 
displacement and stress 
field form. For this 4- 
element mixed-method run, 
the maximum average 
normal stress on the x face 
was found to occur at node 
1 , and was found to be 
approximately 20.3 k.s.i.. 

The average shear stresses 
were all found to be less 
than 1 k.s.i. 

3rd run. 

The third and final run that was performed in NESSUS had a much finer mesh or node 
scheme as is shown in 
Figure 4. Not all elements 
are shown. The beam was 
divided into 40 finite 
elements, the node 
numbers are shown in ()’s. 

There are 20 element edges 
that are equally spaced 
across the beam’s length, 
and 2 element edges that 

are equally spaced across the beam’s height. The element aspect ratio for this beam is 



Figure 4. Global node scheme for 3rd NESSUS run: 
40 element, displacement method. 


Force 

Vector 



[Figure 3. Global node scheme for NESSUS 2nd run: 
4 element, mixed method. 
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1.5:1. This number is the ratio of the element’s maximum dimension to its minimum 
dimension. A conservative approach for any finite element mesh is to try to get this ratio 
to be as close to 1:1 as possible. This model had a total of 63 nodes that would be 
reported on in the output file of NESSUS. For this run, NESSUS used the displacement 
method of analysis to calculate the results. The maximum average normal stress on the x 
face of all of the 63 nodal points within the body was found to be 19.5 k.s.i. This stress 
state was found to occur at node 1. Again, all shear stresses on any face at any of the 
nodal points within the body were found the be less than l k.s.i. 

How do we know that our values are correct? One way to build confidence in the 
answers to any black box is to do simple calculations that give an approximate, quick, and 
cheap answer. The results for the three runs mentioned were checked by using simple 
beam theory. 

Simple Beam Theory 

This theory of analysis is also called technical or Euler’s beam theory. It is based 
on certain assumptions that simplify the deformation pattern so that the strain field for a 
cross section of the member can be determined. The key assumption of this theory is 
that plane sections before deformation (or loading) must remain plane after deformation 
(or loading). This is an exact assumption for axially loaded prismatic bars, circular 
prismatic torsion rods, and for prismatic beams in pure bending. This assumption 
approaches reality when the shear deformation of the beam is negligible. It has been 
shown that the bending deflection is about 100 times more than the shear deflection when 
the length to height ratio of the beam is 10:1; and, as the length to height ratio increases 
the significance of the bending deflection relative to the shear deflection also increases. 
Thus, for the beam to be analyzed, simple beam theory is a justifiable method [3], 

Also, note that there are two loads whose separate effects on the beam must be 
superimposed in order to obtain the effects of both loads acting simultaneously. In order 
to do this, one must have a linear elastic material which remains in the linear elastic range 
after loading; that is, the stresses within the beam must not exceed their yielding value for 
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that material. Also, the beam must exhibit small deformations. These deformations are 
considered small when they produce small strains compared to unity; thus, higher order 
terms involving theses small strains are also neglected [3]. 

The first stress calculated was the normal stress in the x direction due to the axial 
load. This stress was considered constant throughout the beams length and was 
calculated using the equation 


<7 


XX 



Eq. 1 


where is the stress on the x face in the x direction (hence the subscript), P x is the axial 
load shown in Figure 1 , and A is the cross sectional area of the beam. 

The second stress calculated was the bending stress due to the bending moment 
within the body. Note from Figure 1 that there are no pure moments applied to the beam, 
but there does exist a moment at any cross section within the beam, due to the vertical 
load P y and the fact that one end is fixed. The bending stress varies with x and y. It can 

be calculated from the flexure formula: 

M{x)y' 


<T„ = — 


Eq. 2 


where M(x) is the bending moment within a cross section of the beam at x, I is the area 
moment of inertia of the beam, and y' is not the coordinate used in Figure 1 but one that 
has its zero at the neutral axis of the beam. Its positive direction must be the same as that 
of the coordinate y; also, the moment is positive if it produces a positive curvature of the 
beam (smiley face) with respect to the x coordinate shown in Figure 1 and the y' 
coordinate just mentioned. 

When the moments were summed about an arbitrary centroid of a cross section of 
the beam, the following equation resulted 
M(x) = P y L — P y x Eq. 3 

where the only undefined term is L, the length of the beam. Note from equation 2 that the 
maximum stress due to the bending moment occurs when y' is a maximum (at the outer 
surface), and when the bending moment is a maximum (at the fixed end, x=0; as verified 
from equation 3). Therefore, since the normal stress due to the axial load is constant 
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throughout, equations 1, 2, and 3, can be combined to give the following equation that 
represents the maximum total normal stress on the x face 

A bh /n 

where h is the height of the beam, and ^/\2 * s area moment of inertia for this beam. 
Note, this state of stress occurs at x=0, y=0 (/ = “%). After substituting into equation 


4 the numerical values from Figure 1, a value of 20.17 k.s.i. was obtained for the maximum 
total normal stress in the x direction. Recall that the yield strength in tension for this 
type of steel is 36 k.s.i; therefore, the beam remained in the linear elastic range after 
loading. The shear stress was found by using the shear formula 


V 

Txy ~ 21 




Eq.5 


where V is the shear force within a cross section of the beam, I is the area moment of 
inertia of the beam, h is the beam height, and y' is the coordinate system used for the 
flexure formula [1]. From equation 5, one can see that the shear stress is a maximum at 
the neutral axis of the beam, where y is a minimum (zero). When we plug in numerical 
values into equation 5, the result for the maximum shear stress is 0.5 k.s.i, which is 
nowhere close to the yield strength in shear, 21 k.s.i. Both the maximum normal (x dir) 
and the shear stresses were below their respective yield stresses. Therefore, from a 
deterministic standpoint, this beam does not fail by permanently deforming. 


Comparing NESSUS to Simple Beam Theory 

The maximum average normal stress within the body was found to be the most 
significant stress. The magnitude of this stress and its nodal location from all three runs 
in NESSUS, as well as the results from our simple beam theory calculations are shown in 
Table 2. This table also shows the percent difference of each NESSUS run with respect 
to the simple beam theory calculations and the yield strength in tension for ASTM-A36. 
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Table 2. Maximum average normal stress for NESSUS runs and simple team theory calculations. 
The percent difference was calculated with respect to the simple bean theory value. 


Maximum Average 

NESSUS 

NESSUS 

NESSUS 

Simple Beam 

Respective 

Stress 

4 element 

4 element 

40 element 

Theory 

yield 


Displacement 

Mixed 

Displacement 


strengths for 


Method 

Method 

Method 


ASTM-A36 

Node 1, normal 

15 k.s.i 

20.3 k.s.i. 

19.5 k.s.i. 

20.2 k.s.i. 

36 k.s.i. 

% Difference 

25.7 % 

0.5 % 

3.5 % 

0 % 

N/A 


Note from Table 2 that when the amount of elements used to model our beam was 
increased from 4 to 40 (both using the displacement method of analysis), the normal 
stress approached the value obtained from simple beam theory. The percent difference 
was reduced from 25.7 % to 3.5 %. The mixed method of analysis converged to 0.5% of 
the simple beam theory value using only 4 elements. It took about 12 seconds of 
computer time, while the displacement method of analysis needed 40 elements and about 
19 seconds of computer time to get within range of the simple beam theory stress values. 
Also note from the table that none of the maximum average stresses for any run exceeded 
its respective yield strength of the material. Now let us compare NESSUS to a finite 
element code that I wrote. However, before we do that we shall discuss the theory of the 
finite element method by way of weighting the residuals. 


Finite Element Method by Weighting the Residuals 

Using the assumptions mentioned above in the NESSUS portion, the 3D elasticity 
equations were simplified to the following governing differential equations of plane stress: 

roi 


[Lf([C][Lp}) = 


0 


Eq. 8 


where L is a differential operator matrix, C is our linear elastic constitutive model for 
plane stress, and the last term is the displacement vector[2]. Equation 8 can be rewritten 
to show all elements of each tensor as follows 
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Eq. 9 
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where E is the modulus of elasticity of the steel beam (30,000 ksi), and v is Poisson’s 
ratio (0.3). 

Let us begin the method of residual weighting by premultiplying Equation 9 by a weight 
function and integrating over the domain of the problem. The result is the following set of 
equations 

|[M'][Z.] r [CltK«}‘ i£1 = = j°J Eq.I0 

where W is our weight function matrix. Once equation 10 was integrated by parts and we 
realize that there are only nodal forces we end up with the following system of equations 


J([LIW'] r ) r [C][Z.]{«}^ = !' r } Eq.ll 

o 

where F is a vector whose elements are the external forces acting on the system at the 
nodes in for each direction[2]. For example, if a system has 63 nodes and everything 
occurs in either the x or y direction (planar), then F would be a vector with 126 rows or 
entries. After equation 1 1 is descritized by applying it to each element and summing all 
of the element equations together, the next step of the finite element method is to 
approximate our displacement field. For this problem, our displacement fields (u and v) 
were approximated by bilinear interpolation from the nodal values. This approximation is 
mathematically realized by the following equation 




X/V,(x,.vH 

»=l 

£/V,(x,;y)v,. 


i=i 


" N , 0 /V 2 0 /V 3 0 A 4 0 

0 N t 0 N 2 0 A 3 0 N a 






Eq. 12 


where u t and v ( are the nodal displacements of the ith node of the element e in the x and y 
directions, respectively. The N’s are the shape functions for the displacement fields. We 
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must also make the weight function more specific by setting it equal to the transpose of 
the shape function matrix 

W = Eq.13 

hence, because of this step, we have used the Galerkin approach of the finite element 
method[2]. When we substitute equations 12 and 13 into equation 11 , and if we define an 
element strain interpolation matrix as 

[fl] = [L][W] Eq.14 


we end up with the following set of equations 








Eq.15 


where 

[K‘]= J[5] r [C][B]^ Eq.16 

n' 

is defined as the element stiffness matrix and is equal to the term to be integrated (shown 
in brackets) in equation 15. The stiffness matrix now needs to be calculated. The matrix 
multiplication, while time consuming, is not difficult. The integration over a two 
dimensional domain is the difficult part of this step. Gauss quadrature, a typical type of 
numerical integration, will be used to integrate the 64 terms of each element’s 8 by 8 
stiffness matrix. In order to do this, we need to change our global coordinate system 
(x,y), to a local coordinate system (£ and 77 ). The £axis values range from -1 to 1 as does 
the 77 axis. Bilinear maping was used to map the x and y coordinates to the £ and 77 
coordinates. The resulting mapping functions, which are the same shape functions for our 
field variable and our weight function, are given by the following equations 

^2(C- 7 7) = ^(i + C- T 7-Cn) 

Eq.17 

WjCC.t?) = -(i + C + 77 + £*?) 

W 4 (£, 77 ) = -J-(l-C + 77 -^ 77 ) 
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Note that in order for the coordinate mapping to be complete, the limits of integration 
must change from -1 to +1 for both integrals, and our integration variables must change in 
the following manner 


dQ e = dx ■ dy = \ J e • d£ • drj 


Eq.18 


where J is the Jacobian matrix[2], The determinant of this matrix is given by the 
following equation 



dx e dy e dx e dy e 

dri drf dC, 


Eq.19 


where x and y are the global coordinates of the element whose stiffness matrix is being 
calculated, given as a function of £ and r\. The resulting equation for the element stiffness 


matrix comes about by substituting equation 19 into 18, and then 18 into 16. 


[*:']= jj[B] r [C][B] 


- 1-1 


dx e dy e dx e dy £ 


l l 


dC,dr\ = J iMi.vWn 


Eq.20 


-l-i 


dri dr\ , 

where the integrand in equation 20, denoted by ^(£, 77 ), is the function used by Gauss 
quadrature. The integration operators will be replaced by the sum of the product of the 
function evaluated at a certain £ and r\ value and the weight values for those Gauss points. 


Therefore, since 2nd order Guass quadrature is sufficient for each coordinate direction [2], 
equation 20 then becomes 

[V] Eq.21 


g=\h=\ 

where 4Ti — = 0.5773502692, £ 2 = *?2 = -0-5773502692 and the weight functions, w’s, are all 

equal to 1.0. Using equation 21, an element stiffness matrix can be obtained for all 
isoparametric quadrilateral elements. Remember we are trying to analyze a cantilevered beam 
and get results comparable to that of NESSUS and simple beam theory. Therefore, our system 
will be exactly like the 40-element beam modeled in NESSUS, is shown again in Figure 5. 
Therefore, all 40 element stiffness matrices can be obtained; and, since all of the elements look 
the same, all the stiffness matrices will look the same. Each element stiffness matrix will be 
multiplied by a nodal displacement vector to equal an external nodal force vector. Now the 
question is, how is the global stiffness matrix assembled? Well, due to our procedure in solving 
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the plane stress equations, the 
nodal displacement vector for each 
element is a list of displacements 
in the x and y direction, for each 
node in a counter clockwise 
manner. Also, the nodal force 
vector for each element is a list of 
external nodal forces in the x and y 
direction, for each node in a 
counter clockwise manner. For 
example, element 1 of Figure 5. 
shown in Figure 6. has nodes 1 . 2. 
23, and 22, associated with it. 
Therefore, the finite element 
system of equations for element 1 
looks something like 

| ' Kjj 1 j|(5 1 1 = | f ] | , which, when we expand and include (for clarity) the displacement vector 

terms above the stiffness matrix to account for its columns, we arrive at the following set of 
equations 





This can be done for all 40 elements, except the superscript 1 changes to the element 
number, and the subscripts 1, 2, 23, and 22 of the nodal displacement and force vector 
terms respectively change to the 1st, 2nd, 3rd, and 4th node number of the associated 


element. The 1st node of an element is the lower left node. Proceeding in a counter- 
clockwise manner one can then obtain the 2nd through the 4th element node number. 
Note that in equation set 22, any term in the stiffness matrix can be described by its row 
and column, which is in terms of the nodal displacement vector components. For 
example, in equation 22, 

45 = K v2m 2 i Eq.23 

Therefore, when this term is moved to the global set of equations it is placed in the slot 
given by the terms of the right hand side of equation 23. After all the terms of the 


stiffness matrix in equation 22 were assembled in the global set of equations, the result is 


the following incomplete set of equations 


“l 

V, 

u 2 

v 2 • • 

u 22 

v 22 

W 23 

V 23 ' * W 63 V 63 



V ii 

K [ 12 

K l 13 

K' 14 ■ 

■ K l 17 

K' 18 

K \ 5 

AT'tfi ‘ 

U\ 


f/lxl 

K l 21 

K 1 22 

K 1 23 

K'24 

K l 27 

AT'28 

K 1 25 

K' 26 

v l 


fir 

tf l 3J 

K X 32 

K X 33 

K'u 

A' 1 37 

K'u 

**35 

K' 36 

u 2 


fu 

K' 41 

K X 42 

K 5 43 

K' 44 

AT 1 47 

K' 48 


AT'46 

V 2 


fly 

K X 1\ 

K ] 12 

K'j 3 

K' 74 

K'm 

A" 1 78 

K'l 5 

AT '76 

“22 


fllx 

, 

2 

K'n 

K't 4 

K 1 87 

K' 88 

AT *85 

A" 1 86 

v 22 


flly 

tf'si 

K [ 52 

K' 53 

K ] 54 

K ] sn 

K'i 8 

K 1 55 

K'i6 

«23 


fl3x 

K ] 61 

K' 6 2 

K'i 3 

K'w 

K X 67 

K \ 8 

K' 65 

K l 66 

v 23 


fl3y 









«63 


fttx 









. V 63, 


[ ftty 


Eq.24 


where equation 24 becomes the complete set of equations only after this is done in the 


same manner for all 40 elements. If more than one term appears in the same row and 
column slot, then the algebraic sum of the numbers is calculated. 


Once the complete (global) set of equations are assembled, all of the terms of 
equation 24 will either be zero or have a numerical value. The next step would be to 
apply the boundary conditions to the now algebraic differential equation. Recall that for a 
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cantilevered beam the displacements in the x and y directions are zero at the fixed end. 
Referring to Figure 5 one can tell that this implies that the terms w, v, , u 22 , v 22 u 43 v 43 are 

all zero in equation 24. The rest of the terms in the displacement vector { 5 } are the 
unknown displacements of each node in each direction that need to be solved for. Also, 
the force vector {/} , which is the list of external nodal forces for each direction, contains 
unknown terms at the same location where the displacements were known. In other 
words, f\ x ,f \ y , f 2 2x > fiiy - /43X > an d f 43y are the external forces to be calculated. The 

known nodal forces are f 42x and f 42y which are 5001bs and lOOOlbs, respectively. All of 


the other terms in the force vector are zero, this is verified by referring again to Figure 5 
(i.e. the external nodal force in the x direction at node 1 1 is zero — there is no force applied 
there). 


The next step is to reduce the set of equations. This is done by removing the 
rows of the force vector, the stiffness matrix, and the displacement vector that correspond 
to the rows of the displacement vector that have known values (this problem, the ones 
equal to zero). Then the columns of the stiffness matrix that correspond to the known 
terms of the displacement are also removed. This forms our reduced set of equations. 
Using this set the unknown displacements are then calculated in the following manner 


{8 REDUCED } = {/ REDUCED }[ % REDUCED ] Eq.25 

The complete set of nodal displacements (global displacement vector) was then obtained 
by simply expanding its reduced vector in the places where it was reduced and plugging in 
the terms that were taken out (this case all 0’s at nodes 1,22,43). The global force vector 
can then be calculated using the complete set of equations (equation 24, complete) 

[*]{*}={/} 


where K is the global stiffness matrix (before reduction). Note that the only new terms 
that will be obtained are the reaction forces at the fixed end of the beam. 

Once the nodal displacements and forces are known, other needed items like the 
stresses or strains at the nodes can be calculated. For this report, we are concerned with 
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the stresses within the beam. Recall that the first equation we dealt with in this FEM by 
weighted residuals section was the governing equation for plane stress 

[L] r ([ClL]{5}) = j° j 

where L, C and S were expanded in equation 9. This equation comes about by 
substituting the strain-displacement equations of plane stress into the stress-strain 
equations of plane stress and finally that result into the equilibrium equations for plane 
stress. However, we can obtain an equation for the stress as a function of displacement 
by substituting the strain-displacement equations into the stress-strain equations for 
plane stress. The result is the following set of equations 


rr 

E 

°yy 

1-v 2 

x *y. 



1 v 0 
v 1 0 

0 0 izr 
2 . 


I ° f 

d_ | u(x,y) 
dy \v(x,y) 

d d 


Eq.26 


vL dy <^J J 

The displacement vector is then replaced by equation 12 and the result is the following 
equation that can be used to calculate the stress at any point within the beam. 
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Eq.27 


IL^v J 

For the beam shown in Figure 5, the maximum average normal stress in the x direction was 
calculated to be 18.8 k.s.i. It was found to occur at node 1 of the beam. Since the 
maximum stress did not exceed the tensional yield stress of 36 k.s.i. for the material, the 
beam did not fail by permanently deforming. 


Comparing NESSUS to a Finite Element Code 

Since the maximum average normal stress in the x direction was the most 
significant stress within the beam, we shall compare its values obtained from the three 
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NESSUS runs to the value obtained from my finite element code. Figure 7 shows the 
maximum average stress in the x direction obtained from all three runs as well as my finite 
element code. 

Figure7. NESSUS results compared to a deterministic finite element code. 
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The percent difference between the NESSUS 40 element displacement run and the 
deterministic FE code 40 element run was calculated to be 3.59%. Now that confidence 
has been gained with the results from the NESSUS finite element code, let us talk about 
the most important aspect in NESSUS. ..its probabilistic capabilities. 

Probabilistic NESSUS 

NESSUS is a probabilistic finite element code. All of the values that we have 
discussed thus so far have been calculated from the average values of the design variables. 
If you will recall some of the variables had only one deterministic value, but the 
probabilistic design variables had a range of values and associated probabilities of 
occurring. One form of output that a user might request is the probability density 
function (p.d.f.) for the stresses that were calculated. However, for structural 
applications, many times it is desirable that the stress at every point within a material 
stay below that material’s yield stress. When stresses exceed their yield values, plastic or 
permanent deformation has occurred, or will soon occur [3]. Therefore, since many times 
an analyst is only concerned with the chances of the solution (or output, or field) variable 
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being below a certain value, i.e., the stress within the body being less than the material s 
yield strength, NESSUS can be used to calculate a cumulative distribution function, a 
c.d.f., for the requested output. A cumulative distribution function for any random 
variable (this case our stress solution is a r.v. due to the r.v. inputs) is still the same set of 
values, but there is another value associated with the c.d.f. for a random variable. This 
value is the probability or chance that the random variable will be equal to or less than 
that value. The cumulative distribution function for the normal stress in the x direction at 
node 1 for the 4 element mixed method run was obtained from NESSUS and graphed in 
Figure 8. 

Figure 8. Cumulative distribution function for node 1 O ^ mixed method run. 


CUMULATIVE DISTRIBUTION FUNCTION: 
Normal Stress, X Direction, Node 1 



Normal Stress, X Direction, Node 1 (k.s.i.) 


From Figure 8, one can deduce that, due to the uncertainties in our model, there is a 
99.99% probability that the normal, x-directional stress at node 1 will be equal to or less 
than 30.5 k.s.i. Recall that the yield strength in tension for this material is 36 k.s.i. 
Therefore, there is less than a 1/10000 probability that the cantilevered beam analyzed 
will fail by permanently deforming or yielding. 

Conclusion 

The cantilevered beam was analyzed in the probabilistic finite element code NESSUS with 
three different runs. For the 4 element displacement method run, the maximum average 


NASA/CR— 2001-211112 


151 



normal x-direction stress within the beam was found to be 15 k.s.i. For the 4 element 
mixed method run, the maximum average normal x-direction stress was found to be 20.3 
k.s.i. The 40 element displacement method calculated the maximum average normal x- 
direction stress to be 19.5 k.s.i. All maximum normal stresses were below the yield 
strength in tension for the material, which was found to be 36 k.s.i. Also, all maximum 
shear stresses were less than 1 k.s.i, which is well below the material yield strength in 
shear of 25 k.s.i. It was noted that the displacement method needed 40 elements to 
converge upon the value obtained using simple beam theory and while the mixed method 
of analysis needed only 4 elements. Because the normal stress was found to be the most 
significant stress it was calculated using simple beam theory to be 20.2 k.s.i. A 
deterministic finite element code was written and used the displacement method of 
analysis on a 40 element beam. The maximum average normal x-direction stress was 
calculated to be 18.8 k.s.i. The percent difference between the 40-element displacement 
NESSUS run and the other finite element code run was found to be 3.59%. NESSUS was 
observed to yield good results when compared to simple beam theory and another finite 
element code. The probabilistic aspect of NESSUS let us obtain the cumulative 
distribution function for the normal stress in the x direction at the same location of the 
maximum average stress (at the fixed end). From this c.d.f., we observed that the stress at 
this node would be equal to or less than 30.5 k.s.i. 99.99% of the time. Knowing that the 
yield stress in tension for the material was 36 k.s.i., it was concluded that the beam had 
less than a 1/10000 probability of failure by yielding. 
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